Posted By


Bangonkali on 03/24/12

Tagged


Statistics


Viewed 607 times
Favorited by 0 user(s)

Related snippets


Gaussian Elimination


/ Published in: C++
Save to your folder(s)

Numerical Methods application for solving system of equation using Gaussian Elimination based on this Wikipedia article: http://j.mp/GV3PcN


Copy this code and paste it in your HTML
  1. /*
  2.  * File: main.cpp
  3.  * Author: Bangonkali
  4.  *
  5.  * Created on March 23, 2012, 8:14 PM
  6.  */
  7. #include <iostream>
  8. #include <cmath>
  9. #include <cstdlib>
  10.  
  11.  
  12. // example taken from wikipedia
  13. // http://en.wikipedia.org/wiki/Gaussian_elimination#Example
  14.  
  15. using namespace std;
  16.  
  17. //void printMatrix(
  18. // int N,
  19. // float A[20][21]
  20. //)
  21. //{
  22. // for (int i = 0; i < N; i++)
  23. // {
  24. // for (int j = 0; j <= N; j++)
  25. // {
  26. // cout << A[i][j] << '\t';
  27. // }
  28. //
  29. // cout << endl;
  30. // }
  31. //
  32. // cout << endl;
  33. //}
  34.  
  35. int gauss(
  36. int N, // number of unknowns
  37. float A [20] [21], // coefficients and constants
  38. float r[20]
  39. )
  40. {
  41. // printMatrix(N, A);
  42. float multiplier, divider;
  43. bool fin = false;
  44.  
  45. // forward substitution
  46. for (int m=0; m<=N; m++)
  47. {
  48. for (int i=m+1; i<=N; i++)
  49. {
  50. multiplier = A[i][m]; // current row - root
  51. divider = A[m][m]; // first row - root
  52.  
  53. if (divider == 0)
  54. {
  55. return 1;
  56. }
  57.  
  58. for (int j=m; j<=N; j++)
  59. {
  60. if (i == N)
  61. {
  62. break;
  63. fin = true;
  64. }
  65.  
  66. A[i][j] = (A[m][j] * multiplier / divider) - A[i][j];
  67. // cout << "A[" << i << "][" << j << "] = (A[" << m << "][" << j << "] * " << multiplier << " / " << divider << ") - A[" << i << "][" << j << "]" << endl;
  68. // cout << "Current cell " << i << j << endl;
  69. // cout << "Current set " << m << m << endl;
  70. // printMatrix(N, A);
  71. }
  72.  
  73. for (int j=m; j<=N; j++)
  74. {
  75. A[i-1][j] = A[i-1][j] / divider;
  76. }
  77.  
  78. if (fin)
  79. break;
  80. }
  81. }
  82.  
  83. // printMatrix(N, A);
  84.  
  85. // back substitution
  86. float s = 0;
  87. r[N-1] = A[N-1][N];
  88. int y = 0;
  89. for (int i = N-2; i >= 0; i--)
  90. {
  91. s = 0;
  92. y++;
  93. for (int x = 0; x < N; x++)
  94. {
  95. s = s + (A[i][N-1-x] * r[N-(x+1)]);
  96. // cout << "A[" << i << "][" << N-1-x << "] = " << A[i][N-1-x] << " ";
  97. // cout << " r[" << N-(x+1) << "] = " << r[N-(x+1)] << " " << endl;
  98. // cout << "s = " << s << endl;
  99.  
  100. if (y == x+1) break;
  101. }
  102. r[i] = A[i][N] - s;
  103. // cout << "r[" << i << "] = A[" << i << "][" << N << "] - " << s << " = " << r[i] << endl;
  104. }
  105.  
  106. // for (int i = 0; i < N; i++)
  107. // {
  108. // cout << r[i] << endl;
  109. // }
  110.  
  111.  
  112. }
  113. /*
  114.  *
  115.  */
  116. int main(int argc, char** argv) {
  117. float A[20][21];
  118. float X[20];
  119. bool err;
  120. A[0][0] = 2;
  121. A[0][1] = 1;
  122. A[0][2] = -1;
  123. A[0][3] = 8;
  124.  
  125. A[1][0] = -3;
  126. A[1][1] = -1;
  127. A[1][2] = 2;
  128. A[1][3] = -11;
  129.  
  130. A[2][0] = -2;
  131. A[2][1] = 1;
  132. A[2][2] = 2;
  133. A[2][3] = -3;
  134.  
  135. gauss(3, A, X);
  136. for (int i=0; i<3; i++) cout << X[i] << " ";
  137.  
  138. return 0;
  139. }

URL: http://j.mp/GV3PcN

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.