Skip to content

Commit e9527bb

Browse files
committed
Added Gaussian Elimimnation in C++
1 parent d740d52 commit e9527bb

File tree

2 files changed

+158
-3
lines changed

2 files changed

+158
-3
lines changed
Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
#include <iostream>
2+
#include <algorithm>
3+
#include <stdexcept>
4+
#include <vector>
5+
using namespace std;
6+
7+
class EquationSolver
8+
{
9+
private:
10+
11+
// subtract row from another row after multiplcation
12+
void subRow(vector<long double> &toBeSubtracted, const vector<long double> &toSubtract,
13+
long double mult, int start)
14+
{
15+
for (int i = start; i < toBeSubtracted.size(); i++)
16+
toBeSubtracted[i] -= mult * toSubtract[i];
17+
}
18+
// A start variable is given since if we know all values are 0, we don't need
19+
// to subtract.
20+
21+
22+
// swap two rows
23+
void swapRow(vector<long double> &eqn1, vector<long double> &eqn2, int start)
24+
{
25+
for (int i = start; i < eqn1.size(); i++)
26+
swap(eqn1[i], eqn2[i]);
27+
}
28+
// A start variable is given since if we know all values are 0, we don't need
29+
// to swap.
30+
31+
32+
void gaussianElimination(vector<vector<long double>> &eqns, int varc)
33+
{
34+
// 'eqns' is the matrix, 'varc' is no. of vars
35+
for (int i = 0; i < varc - 1; i++)
36+
{
37+
if (eqns[i][i] == 0)
38+
{
39+
for (int j = i + 1; j < varc; j++)
40+
{
41+
if (eqns[j][i] != 0)
42+
{
43+
swapRow(eqns[j], eqns[i], i);
44+
break;
45+
}
46+
}
47+
48+
if (eqns[i][i] == 0)
49+
throw runtime_error("Singular matrix");
50+
// Cannot solve since coefficient of one variable is always zero
51+
}
52+
53+
for (int j = i + 1; j < varc; j++)
54+
if (eqns[j][i])
55+
subRow(eqns[j], eqns[i], (eqns[j][i] / eqns[i][i]), i);
56+
}
57+
58+
if (eqns[varc - 1][varc - 1] == 0)
59+
{
60+
if (eqns[varc - 1][varc] == 0)
61+
throw runtime_error("Singular matrix");
62+
// Cannot solve since final equation is '0*xn = 0'
63+
else
64+
throw runtime_error("No solutions");
65+
// Cannot solve since final equation is '0*xn = c'
66+
}
67+
}
68+
69+
70+
void gaussJordan(vector<vector<long double>> &eqns, int varc)
71+
{
72+
// 'eqns' is the (Row-echelon) matrix, 'varc' is no. of vars
73+
for (int i = varc - 1; i >= 0; i--)
74+
{
75+
eqns[i][varc] /= eqns[i][i];
76+
eqns[i][i] = 1;
77+
for (int j = i - 1; j >= 0; j--)
78+
{
79+
eqns[j][varc] -= eqns[j][i] * eqns[i][varc];
80+
eqns[j][i] = 0;
81+
}
82+
}
83+
}
84+
85+
86+
vector<long double> backSubs(vector<vector<long double>> &eqns, int varc)
87+
{
88+
// 'eqns' is matrix, 'varc' is no. of variables
89+
vector<long double> ans(varc);
90+
for (int i = varc - 1; i >= 0; i--)
91+
{
92+
long double sum = 0;
93+
for (int j = i + 1; j < varc; j++)
94+
sum += eqns[i][j] * ans[j];
95+
96+
ans[i] = (eqns[i][varc] - sum) / eqns[i][i];
97+
}
98+
return ans;
99+
}
100+
101+
public:
102+
vector<long double> solveByGaussJordan(vector< vector<long double> > eqns, int varc)
103+
{
104+
this->gaussianElimination(eqns, varc);
105+
this->gaussJordan(eqns, varc);
106+
107+
vector<long double> ans(varc);
108+
for (int i = 0; i < varc; i++)
109+
ans[i] = eqns[i][varc];
110+
return ans;
111+
}
112+
113+
vector<long double> solveByBacksubs(vector< vector<long double> > eqns, int varc)
114+
{
115+
this->gaussianElimination(eqns, varc);
116+
vector<long double> ans = this->backSubs(eqns, varc);
117+
return ans;
118+
}
119+
};
120+
121+
int main() {
122+
int varc = 3;
123+
vector< vector<long double> > equations { { 2 , 3, 4, 6 },
124+
{ 1 , 2, 3, 4 },
125+
{ 3 ,-4, 0, 10 }};
126+
EquationSolver solver;
127+
vector<long double> ans;
128+
try
129+
{
130+
ans = solver.solveByGaussJordan(equations, varc);
131+
}
132+
catch(runtime_error e)
133+
{
134+
cout << "Error found: " << e.what() << endl;
135+
return -1;
136+
}
137+
138+
cout << "The solution is," << endl
139+
<< "x == " << ans[0] << endl
140+
<< "y == " << ans[1] << endl
141+
<< "z == " << ans[2] << endl;
142+
}

contents/gaussian_elimination/gaussian_elimination.md

Lines changed: 16 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -232,7 +232,7 @@ $$
232232
$$
233233

234234
There are plenty of different strategies you could use to do this, and no one strategy is better than the rest.
235-
One method is to subtract a multiple of the top row from subsequent rows below it such that all values beneath the pivot value are zero.
235+
One method is to subtract a multiple of the top row from subsequent rows below it such that all values beneath the pivot value are zero.
236236
This process might be easier if you swap some rows around first and can be performed for each pivot.
237237

238238
After you get a row echelon matrix, the next step is to find the reduced row echelon form. In other words, we do the following:
@@ -314,6 +314,9 @@ In code, this process might look like this:
314314
{% sample lang="c" %}
315315
[import:5-13, lang:"c_cpp"](code/c/gaussian_elimination.c)
316316
[import:19-34, lang:"c_cpp"](code/c/gaussian_elimination.c)
317+
{% sample lang="cpp" %}
318+
[import:22-29, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
319+
[import:37-51, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
317320
{% sample lang="hs" %}
318321
[import:10-17, lang:"haskell"](code/haskell/gaussianElimination.hs)
319322
[import:44-46, lang:"haskell"](code/haskell/gaussianElimination.hs)
@@ -367,7 +370,7 @@ $$
367370
\left[
368371
\begin{array}{ccc|c}
369372
3 & -4 & 0 & 10 \\
370-
0 & \mathbf{\frac{10}{3}} & \mathbf{3} & \mathbf{\frac{2}{3}}
373+
0 & \mathbf{\frac{10}{3}} & \mathbf{3} & \mathbf{\frac{2}{3}}
371374
\\
372375
2 & 3 & 4 & 6
373376
\end{array}
@@ -384,6 +387,8 @@ Here is what it might look like in code:
384387
[import:32-40, lang:"java"](code/java/GaussianElimination.java)
385388
{% sample lang="c" %}
386389
[import:36-41, lang:"c_cpp"](code/c/gaussian_elimination.c)
390+
{% sample lang="cpp" %}
391+
[import:11-19, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
387392
{% sample lang="hs" %}
388393
[import:19-33, lang:"haskell"](code/haskell/gaussianElimination.hs)
389394
[import:42-42, lang:"haskell"](code/haskell/gaussianElimination.hs)
@@ -403,6 +408,8 @@ When we put everything together, it looks like this:
403408
[import:1-45, lang:"julia"](code/julia/gaussian_elimination.jl)
404409
{% sample lang="c" %}
405410
[import:15-48, lang:"c_cpp"](code/c/gaussian_elimination.c)
411+
{% sample lang="cpp" %}
412+
[import:11-67, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
406413
{% sample lang="rs" %}
407414
[import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs)
408415
{% sample lang="hs" %}
@@ -440,6 +447,8 @@ Here it is in code:
440447
[import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
441448
{% sample lang="c" %}
442449
[import:64-82, lang:"c_cpp"](code/c/gaussian_elimination.c)
450+
{% sample lang="cpp" %}
451+
[import:70-83, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
443452
{% sample lang="rs" %}
444453
This code does not exist yet in rust, so here's Julia code (sorry for the inconvenience)
445454
[import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
@@ -481,6 +490,8 @@ In code, it looks like this:
481490
[import:47-64, lang:"julia"](code/julia/gaussian_elimination.jl)
482491
{% sample lang="c" %}
483492
[import:50-62, lang:"c_cpp"](code/c/gaussian_elimination.c)
493+
{% sample lang="cpp" %}
494+
[import:86-99, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
484495
{% sample lang="rs" %}
485496
[import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs)
486497
{% sample lang="hs" %}
@@ -495,7 +506,7 @@ In code, it looks like this:
495506

496507
## Visual Representation
497508

498-
We have thus far used Gaussian elimination as a method to solve a system of equations; however, there is often a much easier way to find a similar solution simply by plotting each row in our matrix.
509+
We have thus far used Gaussian elimination as a method to solve a system of equations; however, there is often a much easier way to find a similar solution simply by plotting each row in our matrix.
499510
For the case of 2 equations and 2 unknowns, we would plot the two lines corresponding to each equation and the $$(x, y)$$ location of their point of intersection would be the solution for $$x$$ and $$y$$.
500511
Similarly, for the case of 3 equations and 3 unknowns, we would plot 3 planes and the $$(x, y, z)$$ location of their point of intersection would be the solution for $$x$$, $$y$$, and $$z$$.
501512

@@ -536,6 +547,8 @@ There are also plenty of other solvers that do similar things that we will get t
536547
[import, lang:"julia"](code/julia/gaussian_elimination.jl)
537548
{% sample lang="c" %}
538549
[import, lang:"c_cpp"](code/c/gaussian_elimination.c)
550+
{% sample lang="cpp" %}
551+
[import, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
539552
{% sample lang="rs" %}
540553
[import, lang:"rust"](code/rust/gaussian_elimination.rs)
541554
{% sample lang="hs" %}

0 commit comments

Comments
 (0)