71 lines
1.8 KiB
Text
71 lines
1.8 KiB
Text
Solve: procedure options (main); /* 11 January 2014 */
|
|
|
|
declare n fixed binary;
|
|
put ('Program to solve n simultaneous equations of the form Ax = b. Please type n:' );
|
|
get (n);
|
|
|
|
begin;
|
|
declare (A(n, n), b(n), x(n)) float(18);
|
|
declare (SA(n,n), Sb(n)) float (18);
|
|
declare i fixed binary;
|
|
|
|
put skip list ('Please type A:');
|
|
get (a);
|
|
put skip list ('Please type the right-hand sides, b:');
|
|
get (b);
|
|
|
|
SA = A; Sb = b;
|
|
|
|
put skip list ('The equations are:');
|
|
do i = 1 to n;
|
|
put skip edit (A(i,*), b(i)) (f(5), x(1));
|
|
end;
|
|
|
|
call Gauss_elimination (A, b);
|
|
|
|
call Backward_substitution (A, b, x);
|
|
|
|
put skip list ('Solutions:'); put skip data (x);
|
|
|
|
/* Check solutions: */
|
|
put skip list ('Residuals:');
|
|
do i = 1 to n;
|
|
put skip list (sum(SA(i,*) * x(*)) - Sb(i));
|
|
end;
|
|
end;
|
|
|
|
Gauss_elimination: procedure (A, b) options (reorder); /* Triangularise */
|
|
declare (A(*,*), b(*)) float(18);
|
|
declare n fixed binary initial (hbound(A, 1));
|
|
declare (i, j, k) fixed binary;
|
|
declare t float(18);
|
|
|
|
do j = 1 to n;
|
|
do i = j+1 to n; /* For each of the rows beneath the current (pivot) row. */
|
|
t = A(j,j) / A(i,j);
|
|
do k = j+1 to n; /* Subtract a multiple of row i from row j. */
|
|
A(i,k) = A(j,k) - t*A(i,k);
|
|
end;
|
|
b(i) = b(j) - t*b(i); /* ... and the right-hand side. */
|
|
end;
|
|
end;
|
|
end Gauss_elimination;
|
|
|
|
Backward_substitution: procedure (A, b, x) options (reorder);
|
|
declare (A(*,*), b(*), x(*)) float(18);
|
|
declare t float(18);
|
|
declare n fixed binary initial (hbound(A, 1));
|
|
declare (i, j) fixed binary;
|
|
|
|
x(n) = b(n) / a(n,n);
|
|
|
|
do j = n-1 to 1 by -1;
|
|
t = 0;
|
|
do i = j+1 to n;
|
|
t = t + a(j,i)*x(i);
|
|
end;
|
|
x(j) = (b(j) - t) / a(j,j);
|
|
end;
|
|
end Backward_substitution;
|
|
|
|
end Solve;
|