Gaussian elimination
Problem: Solve Ax=b using Gaussian elimination then backwards substitution. A being an n by n matrix. Also, x and b are n by 1 vectors. To improve accuracy, please use partial pivoting and scaling.
ALGOL 68
File: prelude_mat_lib.a68<lang algol68># -*- coding: utf-8 -*- #
MODE SCAL = REAL,
VEC = [0]SCAL, MAT = [0,0]SCAL;
INT lwb vec := 1, upb vec := 0; FORMAT scal repr = $g(-0,real width OVER 2)$; FORMAT sub = $"("$, sep=$", "$, bus = $")"$; FORMAT vec repr = $f(sub)n(upb vec - lwb vec)(f(scal repr)f(sep))f(scal repr)f(bus)$;
- OPerators to swap the values of two VECtors #
PRIO =:= = 1; OP =:= = (REF VEC u, v)VOID:
FOR i TO UPB u DO SCAL scal=u[i]; u[i]:=v[i]; v[i]:=scal OD;
OP +:= = (REF VEC lhs, VEC rhs)REF VEC: (
FOR i TO UPB lhs DO lhs[i] +:= rhs[i] OD; lhs
);
OP -:= = (REF VEC lhs, VEC rhs)REF VEC: (
FOR i TO UPB lhs DO lhs[i] -:= rhs[i] OD; lhs
);
OP *:= = (REF VEC lhs, SCAL rhs)REF VEC: (
FOR i TO UPB lhs DO lhs[i] *:= rhs OD; lhs
);
OP /:= = (REF VEC lhs, SCAL rhs)REF VEC: (
SCAL inv = 1 / rhs; FOR i TO UPB lhs DO lhs[i] *:= inv OD; lhs
);
SKIP</lang>File: prelude_gaussian_elimination.a68<lang algol68># -*- coding: utf-8 -*- #
PROC in situ gaussian elimination = (REF MAT a, REF MAT b)REF MAT: (
FOR row TO UPB a-1 DO INT pivot row := row; SCAL max mult := ABS a[row,row]; FOR sub row FROM row + 1 TO UPB a DO SCAL abs a ij = ABS a[sub row,row]; IF abs a ij>=max mult THEN pivot row := sub row; max mult := abs a ij FI OD; # now we have the best row to pivot, do the actual pivot, both these are pointers # IF row NE pivot row THEN a[pivot row,] =:= a[row,]; # swap/pivot the rows # b[pivot row,] =:= b[row,] FI;
IF ABS a[row,row] < small real THEN raise value error("singular matrix") FI; FOR sub row FROM row+1 TO UPB a DO SCAL mult = a[sub row,row]/a[row,row]; # a[sub row,] -:= mult*a[row,] XXX: "unoptimised" # #DB# a[sub row,row+1:] -:= mult*a[row,row+1:];# XXX: "optimised" # b[sub row,] -:= mult*b[row,] OD OD;
- We have a triangular matrix, at this point we can traverse backwards
up the diagonal calculating x,y&z Converting it initial to a diagonal matrix, then to the identity. #
FOR row FROM UPB a BY -1 TO 1+LWB a DO FOR sub row TO row-1 DO SCAL mult = a[sub row,row]/a[row,row]; IF ABS a[row,row] < small real THEN raise value error("Zero pivot encountered") FI; # a[sub row,row] -:= mult*a[row,row]; XXX: Optimised, this line can be TOTALLY removed # b[sub row,] -:= mult*b[row,] OD OD;
- Now we have a diagonal matrix we can simply divide b
by the values along the diagonal of A. #
FOR sub row TO UPB a DO b[sub row,] /:= a[sub row,sub row] OD; b
);
PROC gaussian elimination = (MAT in a, MAT in b)MAT: (
[UPB in a, 2 UPB in a]SCAL a := in a; [UPB in b, 2 UPB in b]SCAL b := in b; in situ gaussian elimination(a,b)
);
SKIP</lang>File: prelude_exception.a68<lang algol68># -*- coding: utf-8 -*- #
MODE FIXABLE = BOOL; FIXABLE fixed exception = TRUE;
MODE SIMPLEOUTC = [0]UNION(CHAR, STRING, INT, REAL, BOOL, BITS); MODE SIMPLEOUTB = [0]SIMPLEOUTC; MODE SIMPLEOUTA = [0]SIMPLEOUTB; MODE SIMPLEOUT = [0]SIMPLEOUTA;
PROC raise = (SIMPLEOUT message)FIXABLE: ( putf(stand error, ($"Exception:"$, $xg$, message, $l$)); stop );
PROC raise value error = (SIMPLEOUT message)FIXABLE: IF raise(message) NE fixed exception THEN exception value error; FALSE FI;
SKIP</lang>File: postlude_exception.a68<lang algol68>SKIP EXIT
exception too many iterations: exception value error:
stop</lang>File: test_Gaussian_elimination.a68<lang algol68>#!/usr/bin/algol68g-full --script #
- -*- coding: utf-8 -*- #
PR READ "prelude_mat_lib.a68" PR; PR READ "prelude_gaussian_elimination.a68" PR; PR READ "prelude_exception.a68" PR;
[,]REAL 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)
); []REAL b = (-0.01, 0.61, 0.91, 0.99, 0.60, 0.02);
upb vec := UPB b;
printf((vec repr, gaussian elimination(a,TRNSP MAT(b))));
PR READ "postlude_exception.a68" PR</lang>Output:
(-.0100000, 1.6027904, -1.6132031, 1.2454941, -.4909897, .0657607)
C
This modifies A and b in place, which might not be quite desirable. <lang c>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
- define mat_elem(a, y, x, n) (a + ((y) * (n) + (x)))
void swap_row(double *a, double *b, int r1, int r2, int n) { double tmp, *p1, *p2; int i;
if (r1 == r2) return; for (i = 0; i < n; i++) { p1 = mat_elem(a, r1, i, n); p2 = mat_elem(a, r2, i, n); tmp = *p1, *p1 = *p2, *p2 = tmp; } tmp = b[r1], b[r1] = b[r2], b[r2] = tmp; }
void gauss_eliminate(double *a, double *b, double *x, int n) {
- define A(y, x) (*mat_elem(a, y, x, n))
int i, j, col; double max, tmp;
for (col = 0; col < n; col++) { j = col, max = A(j, j);
for (i = col + 1; i < n; i++) if ((tmp = fabs(A(i, col))) > max) j = i, max = tmp;
swap_row(a, b, col, j, n);
for (i = col + 1; i < n; i++) { tmp = A(i, col) / A(col, col); for (j = col + 1; j < n; j++) A(i, j) -= tmp * A(col, j); A(i, col) = 0; b[i] -= tmp * b[col]; } } for (col = n - 1; col >= 0; col--) { tmp = b[col]; for (j = n - 1; j > col; j--) tmp -= x[j] * A(col, j); x[col] = tmp / A(col, col); }
- undef A
}
int main(void) { double 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 }; double b[] = { -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 }; double x[6]; int i;
gauss_eliminate(a, b, x, 6);
for (i = 0; i < 6; i++) printf("%g\n", x[i]);
return 0;
}</lang>
- Output:
-0.01 1.60279 -1.6132 1.24549 -0.49099 0.0657607
J
%. , J's matrix divide verb, directly solves systems of determined and of over-determined linear equations directly. This example J session builds a noisy sine curve on the half circle, fits quintic and quadratic equations, and displays the results of evaluating these polynomials.
<lang J>
f=: 6j2&": NB. formatting verb
sin=: 1&o. NB. verb to evaluate circle function 1, the sine
add_noise=: ] + (* (_0.5 + 0 ?@:#~ #)) NB. AMPLITUDE add_noise SIGNAL
f RADIANS=: o.@:(%~ i.@:>:)5 NB. monadic circle function is pi times 0.00 0.63 1.26 1.88 2.51 3.14
f SINES=: sin RADIANS 0.00 0.59 0.95 0.95 0.59 0.00
f NOISY_SINES=: 0.1 add_noise SINES _0.01 0.61 0.91 0.99 0.60 0.02
A=: (^/ i.@:#) RADIANS NB. A is the quintic coefficient matrix
NB. display the equation to solve (f A) ; 'x' ; '=' ; f@:,. NOISY_SINES
┌────────────────────────────────────┬─┬─┬──────┐ │ 1.00 0.00 0.00 0.00 0.00 0.00│x│=│ _0.01│ │ 1.00 0.63 0.39 0.25 0.16 0.10│ │ │ 0.61│ │ 1.00 1.26 1.58 1.98 2.49 3.13│ │ │ 0.91│ │ 1.00 1.88 3.55 6.70 12.62 23.80│ │ │ 0.99│ │ 1.00 2.51 6.32 15.88 39.90100.28│ │ │ 0.60│ │ 1.00 3.14 9.87 31.01 97.41306.02│ │ │ 0.02│ └────────────────────────────────────┴─┴─┴──────┘
f QUINTIC_COEFFICIENTS=: NOISY_SINES %. A NB. %. solves the linear system _0.01 1.71 _1.88 1.48 _0.58 0.08
quintic=: QUINTIC_COEFFICIENTS&p. NB. verb to evaluate the polynomial
NB. %. also solves the least squares fit for overdetermined system quadratic=: (NOISY_SINES %. (^/ i.@:3:) RADIANS)&p. NB. verb to evaluate quadratic. quadratic
_0.0200630695393961729 1.26066877804926536 _0.398275112136019516&p.
NB. The quintic is agrees with the noisy data, as it should f@:(NOISY_SINES ,. sin ,. quadratic ,. quintic) RADIANS _0.01 0.00 _0.02 _0.01 0.61 0.59 0.61 0.61 0.91 0.95 0.94 0.91 0.99 0.95 0.94 0.99 0.60 0.59 0.63 0.60 0.02 0.00 0.01 0.02
f MID_POINTS=: (+ -:@:(-/@:(2&{.)))RADIANS _0.31 0.31 0.94 1.57 2.20 2.83
f@:(sin ,. quadratic ,. quintic) MID_POINTS _0.31 _0.46 _0.79 0.31 0.34 0.38 0.81 0.81 0.77 1.00 0.98 1.00 0.81 0.83 0.86 0.31 0.36 0.27
</lang>
Mathematica
<lang Mathematica>GaussianElimination[A_?MatrixQ, b_?VectorQ] := Last /@ RowReduce[Flatten /@ Transpose[{A, b}]]</lang>
MATLAB
<lang MATLAB> function [ x ] = GaussElim( A, b)
% Ensures A is n by n sz = size(A); if sz(1)~=sz(2)
fprintf('A is not n by n\n'); clear x; return;
end
n = sz(1);
% Ensures b is n x 1. if n~=sz(1)
fprintf('b is not 1 by n.\n'); return
end
x = zeros(n,1); aug = [A b]; tempmatrix = aug;
for i=2:sz(1)
% Find maximum of row and divide by the maximum tempmatrix(1,:) = tempmatrix(1,:)/max(tempmatrix(1,:)); % Finds the maximum in column temp = find(abs(tempmatrix) - max(abs(tempmatrix(:,1)))); if length(temp)>2 for j=1:length(temp)-1 if j~=temp(j) maxi = j; %maxi = column number of maximum break; end end else % length(temp)==2 maxi=1; end % Row swap if maxi is not 1 if maxi~=1 temp = tempmatrix(maxi,:); tempmatrix(maxi,:) = tempmatrix(1,:); tempmatrix(1,:) = temp; end % Row reducing for j=2:length(tempmatrix)-1 tempmatrix(j,:) = tempmatrix(j,:)-tempmatrix(j,1)/tempmatrix(1,1)*tempmatrix(1,:); if tempmatrix(j,j)==0 || isnan(tempmatrix(j,j)) || abs(tempmatrix(j,j))==Inf fprintf('Error: Matrix is singular.\n'); clear x; return end end aug(i-1:end,i-1:end) = tempmatrix; % Decrease matrix size tempmatrix = tempmatrix(2:end,2:end);
end
% Backwards Substitution x(end) = aug(end,end)/aug(end,end-1); for i=n-1:-1:1
x(i) = (aug(i,end)-dot(aug(i,1:end-1),x))/aug(i,i);
end
end </lang>
PARI/GP
If A and B have floating-point numbers (t_REAL
s) then the following uses Gaussian elimination:
<lang parigp>matsolve(A,B)</lang>
If the entries are integers, then p-adic lifting (Dixon 1982) is used instead.
Tcl
<lang tcl>package require math::linearalgebra
set 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}
} set b {-0.01 0.61 0.91 0.99 0.60 0.02} puts -nonewline [math::linearalgebra::show [math::linearalgebra::solveGauss $A $b] "%.2f"]</lang>
- Output:
-0.01 1.60 -1.61 1.25 -0.49 0.07