43 lines
1.6 KiB
Text
43 lines
1.6 KiB
Text
proc Gauss(A, B, X, N); \Gaussian elimination returns [X] for [A]*[X] = [B]
|
|
real A, B, X; int N; \matrix, vector, output vector, number of rows
|
|
int Diag, MaxRow, Row, Col, I;
|
|
real Max, Temp;
|
|
[for Diag:= 0 to N-1 do \partial pivoting uses largest magnitude coefficient
|
|
[MaxRow:= Diag; \ to improve numerical stability
|
|
Max:= A(Diag, Diag);
|
|
for Row:= Diag+1 to N-1 do
|
|
[Temp:= abs(A(Row, Diag));
|
|
if Temp > Max then
|
|
[MaxRow:= Row; Max:= Temp];
|
|
];
|
|
Temp:= A(Diag); A(Diag):= A(MaxRow); A(MaxRow):= Temp; \swap rows
|
|
Temp:= B(Diag); B(Diag):= B(MaxRow); B(MaxRow):= Temp;
|
|
for Row:= Diag+1 to N-1 do
|
|
[Temp:= A(Row, Diag) / A(Diag, Diag); \divide by pivot element
|
|
for Col:= Diag+1 to N-1 do
|
|
A(Row, Col):= A(Row, Col) - Temp*A(Diag, Col);
|
|
A(Row, Diag):= 0.;
|
|
B(Row):= B(Row) - Temp*B(Diag);
|
|
];
|
|
]; \reduced row echelon form
|
|
for Row:= N-1 downto 0 do \back substitution makes [X]
|
|
[Temp:= B(Row);
|
|
for I:= N-1 downto Row+1 do
|
|
Temp:= Temp - X(I)*A(Row, I);
|
|
X(Row):= Temp / A(Row, Row);
|
|
];
|
|
];
|
|
|
|
real A, B, X(6);
|
|
int I;
|
|
[A:= [ [1.00, 0.00, 0.00, 0.00, 0.00, 0.00],
|
|
[1.00, 0.63, 0.39, 0.25, 0.16, 0.10],
|
|
[1.00, 1.26, 1.58, 1.98, 2.49, 3.13],
|
|
[1.00, 1.88, 3.55, 6.70, 12.62, 23.80],
|
|
[1.00, 2.51, 6.32, 15.88, 39.90, 100.28],
|
|
[1.00, 3.14, 9.87, 31.01, 97.41, 306.02] ];
|
|
B:= [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02];
|
|
Gauss(A, B, X, 6);
|
|
for I:= 0 to 6-1 do
|
|
[RlOut(0, X(I)); CrLf(0)];
|
|
]
|