Gaussian elimination
You are encouraged to solve this task according to the task description, using any language you may know.
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_exception.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT PROVIDES
MODE FIXED; INT fixed exception, unfixed exception; PROC (STRING message) FIXED raise, raise value error
END COMMENT
- Note: ℵ indicates attribute is "private", and
should not be used outside of this prelude #
MODE FIXED = BOOL; # if an exception is detected, can it be fixed "on-site"? # FIXED fixed exception = TRUE, unfixed exception = FALSE;
MODE #ℵ#SIMPLEOUTV = [0]UNION(CHAR, STRING, INT, REAL, BOOL, BITS); MODE #ℵ#SIMPLEOUTM = [0]#ℵ#SIMPLEOUTV; MODE #ℵ#SIMPLEOUTT = [0]#ℵ#SIMPLEOUTM; MODE SIMPLEOUT = [0]#ℵ#SIMPLEOUTT;
PROC raise = (#ℵ#SIMPLEOUT message)FIXED: (
putf(stand error, ($"Exception:"$, $xg$, message, $l$)); stop
);
PROC raise value error = (#ℵ#SIMPLEOUT message)FIXED:
IF raise(message) NE fixed exception THEN exception value error; FALSE FI;
SKIP</lang>File: prelude_mat_lib.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT PRELUDE REQUIRES
MODE SCAL = REAL; FORMAT scal repr = real repr # and various SCAL OPerators #
END COMMENT
COMMENT PRELUDE PROIVIDES
MODE VEC, MAT; OP :=:, -:=, +:=, *:=, /:=; FORMAT sub, sep, bus; FORMAT vec repr, mat repr
END COMMENT
- Note: ℵ indicates attribute is "private", and
should not be used outside of this prelude #
INT #ℵ#lwb vec := 1, #ℵ#upb vec := 0; INT #ℵ#lwb mat := 1, #ℵ#upb mat := 0; MODE VEC = [lwb vec:upb vec]SCAL,
MAT = [lwb mat:upb mat,lwb vec:upb vec]SCAL;
FORMAT sub := $"( "$, sep := $", "$, bus := $")"$, nl:=$lxx$; FORMAT vec repr := $f(sub)n(upb vec - lwb vec)(f(scal repr)f(sep))f(scal repr)f(bus)$; FORMAT mat repr := $f(sub)n(upb mat - lwb mat)(f( vec repr)f(nl))f( vec repr)f(bus)$;
- OPerators to swap the contents 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; # multiplication is faster # FOR i TO UPB lhs DO lhs[i] *:= inv OD; lhs
);
SKIP</lang>File: prelude_gaussian_elimination.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT PRELUDE REQUIRES
MODE SCAL = REAL, REAL near min scal = min real ** 0.99, MODE VEC = []REAL, MODE MAT = [,]REAL, FORMAT scal repr = real repr, and various OPerators of MAT and VEC
END COMMENT
COMMENT PRELUDE PROVIDES
PROC(MAT a, b)MAT gaussian elimination; PROC(REF MAT a, b)REF MAT in situ gaussian elimination
END COMMENT
- using Gaussian elimination, find x where A*x = b #
PROC in situ gaussian elimination = (REF MAT a, b)REF MAT: (
- Note: a and b are modified "in situ", and b is returned as x #
FOR diag TO UPB a-1 DO INT pivot row := diag; SCAL pivot factor := ABS a[diag,diag]; FOR row FROM diag + 1 TO UPB a DO # Full pivoting # SCAL abs a diag = ABS a[row,diag]; IF abs a diag>=pivot factor THEN pivot row := row; pivot factor := abs a diag FI OD; # now we have the "best" diag to full pivot, do the actual pivot # IF diag NE pivot row THEN # a[pivot row,] =:= a[diag,]; XXX: unoptimised # #DB# a[pivot row,diag:] =:= a[diag,diag:]; # XXX: optimised # b[pivot row,] =:= b[diag,] # swap/pivot the diags of a & b # FI;
IF ABS a[diag,diag] <= near min scal THEN raise value error("singular matrix") FI; SCAL a diag reciprocal := 1 / a[diag, diag];
FOR row FROM diag+1 TO UPB a DO SCAL factor = a[row,diag] * a diag reciprocal; # a[row,] -:= factor * a[diag,] XXX: "unoptimised" # #DB# a[row,diag+1:] -:= factor * a[diag,diag+1:];# XXX: "optimised" # b[row,] -:= factor * b[diag,] OD OD;
- We have a triangular matrix, at this point we can traverse backwards
up the diagonal calculating b\A Converting it initial to a diagonal matrix, then to the identity. #
FOR diag FROM UPB a BY -1 TO 1+LWB a DO
IF ABS a[diag,diag] <= near min scal THEN raise value error("Zero pivot encountered?") FI; SCAL a diag reciprocal = 1 / a[diag,diag];
FOR row TO diag-1 DO SCAL factor = a[row,diag] * a diag reciprocal; # a[row,diag] -:= factor * a[diag,diag]; XXX: "unoptimised" so remove # #DB# b[row,] -:= factor * b[diag,] OD;
- Now we have only diagonal elements we can simply divide b
by the values along the diagonal of A. # b[diag,] *:= a diag reciprocal OD;
b # EXIT #
);
PROC gaussian elimination = (MAT in a, in b)MAT: (
- Note: a and b are cloned and not modified "in situ" #
[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: postlude_exception.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT POSTLUDE PROIVIDES
PROC VOID exception too many iterations, exception value error;
END COMMENT
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_exception.a68" PR;
- define the attributes of the scalar field being used #
MODE SCAL = REAL; FORMAT scal repr = $g(-0,real width)$;
- create "near min scal" as is scales better then small real #
SCAL near min scal = min real ** 0.99;
PR READ "prelude_mat_lib.a68" PR; PR READ "prelude_gaussian_elimination.a68" PR;
MAT 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));
VEC b = (-0.01, 0.61, 0.91, 0.99, 0.60, 0.02);
[UPB b,1]SCAL col b; col b[,1]:= b;
upb vec := 2 UPB a;
printf((vec repr, gaussian elimination(a,col b)));
PR READ "postlude_exception.a68" PR</lang>Output:
( -.010000000000002, 1.602790394502130, -1.613203059905640, 1.245494121371510, -.490989719584686, .065760696175236)
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, row, max_row,dia; double max, tmp;
for (dia = 0; dia < n; dia++) { max_row = dia, max = A(dia, dia);
for (row = dia + 1; row < n; row++) if ((tmp = fabs(A(row, dia))) > max) max_row = row, max = tmp;
swap_row(a, b, dia, max_row, n);
for (row = dia + 1; row < n; row++) { tmp = A(row, dia) / A(dia, dia); for (col = dia+1; col < n; col++) A(row, col) -= tmp * A(dia, col); A(row, dia) = 0; b[row] -= tmp * b[dia]; } } for (row = n - 1; row >= 0; row--) { tmp = b[row]; for (j = n - 1; j > row; j--) tmp -= x[j] * A(row, j); x[row] = tmp / A(row, row); }
- 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
Common Lisp
<lang CommonLisp> (defmacro mapcar-1 (fn n list)
"Maps a function of two parameters where the first one is fixed, over a list" `(mapcar #'(lambda (l) (funcall ,fn ,n l)) ,list) )
(defun gauss (m)
(labels ((redc (m) ; Reduce to triangular form (if (null (cdr m)) m (cons (car m) (mapcar-1 #'cons 0 (redc (mapcar #'cdr (mapcar #'(lambda (r) (mapcar #'- (mapcar-1 #'* (caar m) r) (mapcar-1 #'* (car r) (car m)))) (cdr m)))))) )) (rev (m) ; Reverse each row except the last element (reverse (mapcar #'(lambda (r) (append (reverse (butlast r)) (last r))) m)) )) (catch 'result (let ((m1 (redc (rev (redc m))))) (reverse (mapcar #'(lambda (r) (let ((pivot (find-if-not #'zerop r))) (if pivot (/ (car (last r)) pivot) (throw 'result 'singular)))) m1)) ))))
</lang>
- Output:
(setq m1 '((1.00 0.00 0.00 0.00 0.00 0.00 -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.90 100.28 0.60) (1.00 3.14 9.87 31.01 97.41 306.02 0.02) )) (gauss m1) => (-0.009999999 1.6027923 -1.6132091 1.2455008 -0.4909925 0.06576109)
D
<lang d>import std.stdio, std.math, std.algorithm, std.range, std.numeric,
std.typecons;
Tuple!(double[],"x", string,"err") gaussPartial(in double[][] a0, in double[] b0) pure /*nothrow*/ in {
assert(a0.length == a0[0].length); assert(a0.length == b0.length); assert(a0.all!(row => row.length == a0[0].length));
} body {
enum eps = 1e-6; immutable m = b0.length;
// Make augmented matrix. //auto a = a0.zip(b0).map!(c => c[0] ~ c[1]).array; // Not mutable. auto a = a0.zip(b0).map!(c => [] ~ c[0] ~ c[1]).array;
// Wikipedia algorithm from Gaussian elimination page, // produces row-eschelon form. foreach (immutable k; 0 .. a.length) { // Find pivot for column k and swap. a[k .. m].minPos!((x, y) => x[k] > y[k]).front.swap(a[k]); if (a[k][k].abs < eps) return typeof(return)(null, "singular");
// Do for all rows below pivot. foreach (immutable i; k + 1 .. m) { // Do for all remaining elements in current row. a[i][k+1 .. m+1] -= a[k][k+1 .. m+1] * (a[i][k] / a[k][k]);
a[i][k] = 0; // Fill lower triangular matrix with zeros. } }
// End of WP algorithm. Now back substitute to get result. auto x = new double[m]; foreach_reverse (immutable i; 0 .. m) x[i] = (a[i][m] - a[i][i+1 .. m].dotProduct(x[i+1 .. m])) / a[i][i];
return typeof(return)(x, null);
}
void main() {
// The test case result is correct to this tolerance. enum eps = 1e-13;
// Common RC example. Result computed with rational arithmetic // then converted to double, and so should be about as close to // correct as double represention allows. immutable 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]]; immutable b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02];
immutable r = gaussPartial(a, b); if (!r.err.empty) return writeln("Error: ", r.err); r.x.writeln;
immutable result = [-0.01, 1.602790394502114, -1.6132030599055613, 1.2454941213714368, -0.4909897195846576, 0.065760696175232]; foreach (immutable i, immutable xi; result) if (abs(r.x[i] - xi) > eps) return writeln("Out of tolerance: ", r.x[i], " ", xi);
}</lang>
- Output:
[-0.01, 1.60279, -1.6132, 1.24549, -0.49099, 0.0657607]
Delphi
<lang Delphi>program GuassianElimination;
// Modified from: // R. Sureshkumar (10 January 1997) // Gregory J. McRae (22 October 1997) // http://web.mit.edu/10.001/Web/Course_Notes/Gauss_Pivoting.c
{$APPTYPE CONSOLE}
{$R *.res}
uses
System.SysUtils;
type
TMatrix = class private _r, _c : integer; data : array of TDoubleArray; function getValue(rIndex, cIndex : integer): double; procedure setValue(rIndex, cIndex : integer; value: double); public constructor Create (r, c : integer); destructor Destroy; override;
property r : integer read _r; property c : integer read _c; property value[rIndex, cIndex: integer]: double read getValue write setValue; default; end;
constructor TMatrix.Create (r, c : integer);
begin
inherited Create; self.r := r; self.c := c; setLength (data, r, c);
end;
destructor TMatrix.Destroy; begin
data := nil; inherited;
end;
function TMatrix.getValue(rIndex, cIndex: Integer): double; begin
Result := data[rIndex-1, cIndex-1]; // 1-based array
end;
procedure TMatrix.setValue(rIndex, cIndex : integer; value: double); begin
data[rIndex-1, cIndex-1] := value; // 1-based array
end;
// Solve A x = b procedure gauss (A, b, x : TMatrix); var rowx : integer;
i, j, k, n, m : integer; amax, xfac, temp, temp1 : double;
begin
rowx := 0; // Keep count of the row interchanges n := A.r; for k := 1 to n - 1 do begin amax := abs (A[k,k]); m := k; // Find the row with largest pivot for i := k + 1 to n do begin xfac := abs (A[i,k]); if xfac > amax then begin amax := xfac; m := i; end; end;
if m <> k then begin // Row interchanges rowx := rowx+1; temp1 := b[k,1]; b[k,1] := b[m,1]; b[m,1] := temp1; for j := k to n do begin temp := a[k,j]; a[k,j] := a[m,j]; a[m,j] := temp; end; end;
for i := k+1 to n do begin xfac := a[i, k]/a[k, k]; for j := k+1 to n do a[i,j] := a[i,j]-xfac*a[k,j]; b[i,1] := b[i,1] - xfac*b[k,1] end; end;
// Back substitution for j := 1 to n do begin k := n-j + 1; x[k,1] := b[k,1]; for i := k+1 to n do begin x[k,1] := x[k,1] - a[k,i]*x[i,1]; end; x[k,1] := x[k,1]/a[k,k]; end;
end;
var A, b, x : TMatrix;
begin
try // Could have been done with simple arrays rather than a specific TMatrix class A := TMatrix.Create (4,4); // Note ideal but use TMatrix to define the vectors as well b := TMatrix.Create (4,1); x := TMatrix.Create (4,1);
A[1,1] := 2; A[1,2] := 1; A[1,3] := 0; A[1,4] := 0; A[2,1] := 1; A[2,2] := 1; A[2,3] := 1; A[2,4] := 0; A[3,1] := 0; A[3,2] := 1; A[3,3] := 2; A[3,4] := 1; A[4,1] := 0; A[3,2] := 0; A[4,3] := 1; A[4,4] := 2;
b[1,1] := 2; b[2,1] := 1; b[3,1] := 4; b[4,1] := 8;
gauss (A, b, x);
writeln (x[1,1]:5:2); writeln (x[2,1]:5:2); writeln (x[3,1]:5:2); writeln (x[4,1]:5:2);
readln; except on E: Exception do Writeln(E.ClassName, ': ', E.Message); end;
end.
</lang>
- Output:
1.00, 0.00, 0.00, 4.00
Fortran
Gaussian Elimination with partial pivoting using augmented matrix <lang fortran>
program ge
real, allocatable :: a(:,:),b(:) a = reshape( & [1.0, 1.00, 1.00, 1.00, 1.00, 1.00, & 0.0, 0.63, 1.26, 1.88, 2.51, 3.14, & 0.0, 0.39, 1.58, 3.55, 6.32, 9.87, & 0.0, 0.25, 1.98, 6.70, 15.88, 31.01, & 0.0, 0.16, 2.49, 12.62, 39.90, 97.41, & 0.0, 0.10, 3.13, 23.80, 100.28, 306.02], [6,6] ) b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02] print'(f15.7)',solve_wbs(ge_wpp(a,b))
contains function solve_wbs(u) result(x) ! solve with backward substitution real :: u(:,:) integer :: i,n real , allocatable :: x(:) n = size(u,1) allocate(x(n)) forall (i=n:1:-1) x(i) = ( u(i,n+1) - sum(u(i,i+1:n)*x(i+1:n)) ) / u(i,i) end function
function ge_wpp(a,b) result(u) ! gaussian eliminate with partial pivoting real :: a(:,:),b(:),upi integer :: i,j,n,p real , allocatable :: u(:,:) n = size(a,1) u = reshape( [a,b], [n,n+1] ) do j=1,n p = maxloc(abs(u(j:n,j)),1) + j-1 ! maxloc returns indices between (1,n-j+1) if (p /= j) u([p,j],j) = u([j,p],j) u(j+1:,j) = u(j+1:,j)/u(j,j) do i=j+1,n+1 upi = u(p,i) if (p /= j) u([p,j],i) = u([j,p],i) u(j+1:n,i) = u(j+1:n,i) - upi*u(j+1:n,j) end do end do end function
end program
</lang>
Go
Gaussian elimination with partial pivoting by pseudocode on WP page "Gaussian elimination." Scaling not implemented. <lang go>package main
import (
"errors" "fmt" "log" "math"
)
type testCase struct {
a [][]float64 b []float64 x []float64
}
var tc = testCase{
// common RC example. Result x computed with rational arithmetic then // converted to float64, and so should be about as close to correct as // float64 represention allows. a: [][]float64{ {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: []float64{-0.01, 0.61, 0.91, 0.99, 0.60, 0.02}, x: []float64{-0.01, 1.602790394502114, -1.6132030599055613, 1.2454941213714368, -0.4909897195846576, 0.065760696175232},
}
// result from above test case turns out to be correct to this tolerance. const ε = 1e-13
func main() {
x, err := GaussPartial(tc.a, tc.b) if err != nil { log.Fatal(err) } fmt.Println(x) for i, xi := range x { if math.Abs(tc.x[i]-xi) > ε { log.Println("out of tolerance") log.Fatal("expected", tc.x) } }
}
func GaussPartial(a0 [][]float64, b0 []float64) ([]float64, error) {
// make augmented matrix m := len(b0) a := make([][]float64, m) for i, ai := range a0 { row := make([]float64, m+1) copy(row, ai) row[m] = b0[i] a[i] = row } // WP algorithm from Gaussian elimination page // produces row-eschelon form for k := range a { // Find pivot for column k: iMax := k max := math.Abs(a[k][k]) for i := k + 1; i < m; i++ { if abs := math.Abs(a[i][k]); abs > max { iMax = i max = abs } } if a[iMax][k] == 0 { return nil, errors.New("singular") } // swap rows(k, i_max) a[k], a[iMax] = a[iMax], a[k] // Do for all rows below pivot: for i := k + 1; i < m; i++ { // Do for all remaining elements in current row: for j := k + 1; j <= m; j++ { a[i][j] -= a[k][j] * (a[i][k] / a[k][k]) } // Fill lower triangular matrix with zeros: a[i][k] = 0 } } // end of WP algorithm. // now back substitute to get result. x := make([]float64, m) for i := m - 1; i >= 0; i-- { x[i] = a[i][m] for j := i + 1; j < m; j++ { x[i] -= a[i][j] * x[j] } x[i] /= a[i][i] } return x, nil
}</lang>
- Output:
[-0.01 1.6027903945020987 -1.613203059905494 1.245494121371364 -0.49098971958462834 0.06576069617522803]
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>
JavaScript
From Numerical Recipes in C: <lang javascript>// Lower Upper Solver function lusolve(A, b, update) { var lu = ludcmp(A, update) if (lu === undefined) return // Singular Matrix! return lubksb(lu, b, update) }
// Lower Upper Decomposition function ludcmp(A, update) { // A is a matrix that we want to decompose into Lower and Upper matrices. var d = true var n = A.length var idx = new Array(n) // Output vector with row permutations from partial pivoting var vv = new Array(n) // Scaling information
for (var i=0; i<n; i++) { var max = 0 for (var j=0; j<n; j++) { var temp = Math.abs(A[i][j]) if (temp > max) max = temp } if (max == 0) return // Singular Matrix! vv[i] = 1 / max // Scaling }
if (!update) { // make a copy of A var Acpy = new Array(n) for (var i=0; i<n; i++) { var Ai = A[i] Acpyi = new Array(Ai.length) for (j=0; j<Ai.length; j+=1) Acpyi[j] = Ai[j] Acpy[i] = Acpyi } A = Acpy }
var tiny = 1e-20 // in case pivot element is zero for (var i=0; ; i++) { for (var j=0; j<i; j++) { var sum = A[j][i] for (var k=0; k<j; k++) sum -= A[j][k] * A[k][i]; A[j][i] = sum } var jmax = 0 var max = 0; for (var j=i; j<n; j++) { var sum = A[j][i] for (var k=0; k<i; k++) sum -= A[j][k] * A[k][i]; A[j][i] = sum var temp = vv[j] * Math.abs(sum) if (temp >= max) { max = temp jmax = j } } if (i <= jmax) { for (var j=0; j<n; j++) { var temp = A[jmax][j] A[jmax][j] = A[i][j] A[i][j] = temp } d = !d; vv[jmax] = vv[i] } idx[i] = jmax; if (i == n-1) break; var temp = A[i][i] if (temp == 0) A[i][i] = temp = tiny temp = 1 / temp for (var j=i+1; j<n; j++) A[j][i] *= temp } return {A:A, idx:idx, d:d} }
// Lower Upper Back Substitution function lubksb(lu, b, update) { // solves the set of n linear equations A*x = b. // lu is the object containing A, idx and d as determined by the routine ludcmp. var A = lu.A var idx = lu.idx var n = idx.length
if (!update) { // make a copy of b var bcpy = new Array(n) for (var i=0; i<b.length; i+=1) bcpy[i] = b[i] b = bcpy }
for (var ii=-1, i=0; i<n; i++) { var ix = idx[i] var sum = b[ix] b[ix] = b[i] if (ii > -1) for (var j=ii; j<i; j++) sum -= A[i][j] * b[j] else if (sum) ii = i b[i] = sum } for (var i=n-1; i>=0; i--) { var sum = b[i] for (var j=i+1; j<n; j++) sum -= A[i][j] * b[j] b[i] = sum / A[i][i] } return b // solution vector x }
document.write( lusolve( [ [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]
],
[-0.01, 0.61, 0.91, 0.99, 0.60, 0.02]
) )</lang>
- Output:
-0.01000000000000004, 1.6027903945021095, -1.6132030599055475, 1.2454941213714232, -0.4909897195846526, 0.06576069617523138
Julia
Using built-in LAPACK-based linear solver (which employs partial-pivoted Gaussian elimination): <lang julia>x = A \ b</lang>
Mathematica / Wolfram Language
<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>
OCaml
The OCaml stdlib is fairly lean, so these stand-alone solutions often need to include support functions which would be part of a codebase, like these... <lang OCaml> module Array = struct
include Array (* Computes: f a.(0) + f a.(1) + ... where + is 'g'. *) let foldmap g f a = let n = Array.length a in let rec aux acc i = if i >= n then acc else aux (g acc (f a.(i))) (succ i) in aux (f a.(0)) 1
(* like the stdlib fold_left, but also provides index to f *) let foldi_left f x a = let r = ref x in for i = 0 to length a - 1 do r := f i !r (unsafe_get a i) done; !r
end
let foldmap_range g f (a,b) =
let rec aux acc n = let n = succ n in if n > b then acc else aux (g acc (f n)) n in aux (f a) a
let fold_range f init (a,b) =
let rec aux acc n = if n > b then acc else aux (f acc n) (succ n) in aux init a
</lang> The solver: <lang OCaml> (* Some less-general support functions for 'solve'. *) let swap_elem m i j = let x = m.(i) in m.(i) <- m.(j); m.(j) <- x let maxtup a b = if (snd a) > (snd b) then a else b let augmented_matrix m b =
Array.(init (length m) ( fun i -> append m.(i) [|b.(i)|] ))
(* Solve Ax=b for x, using gaussian elimination with scaled partial pivot,
* and then back-substitution of the resulting row-echelon matrix. *)
let solve m b =
let n = Array.length m in let n' = pred n in (* last index = n-1 *) let s = Array.(map (foldmap max abs_float) m) in (* scaling vector *) let a = augmented_matrix m b in
for k = 0 to pred n' do (* Scaled partial pivot, to preserve precision *) let pair i = (i, abs_float a.(i).(k) /. s.(i)) in let i_max,v = foldmap_range maxtup pair (k,n') in if v < epsilon_float then failwith "Matrix is singular."; swap_elem a k i_max; swap_elem s k i_max;
(* Eliminate one column *) for i = succ k to n' do let tmp = a.(i).(k) /. a.(k).(k) in for j = succ k to n do a.(i).(j) <- a.(i).(j) -. tmp *. a.(k).(j); done done done;
(* Backward substitution; 'b' is in the 'nth' column of 'a' *) let x = Array.copy b in (* just a fresh array of the right size and type *) for i = n' downto 0 do let minus_dprod t j = t -. x.(j) *. a.(i).(j) in x.(i) <- fold_range minus_dprod a.(i).(n) (i+1,n') /. a.(i).(i); done; x
</lang> Example data... <lang OCaml> let 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 |] |]
let b = [| -0.01; 0.61; 0.91; 0.99; 0.60; 0.02 |] </lang> In the REPL, the solution is: <lang OCaml>
- let x = solve a b;;
val x : float array = [|-0.0100000000000000991; 1.60279039450210536; -1.61320305990553226;
1.24549412137140547; -0.490989719584644546; 0.0657606961752301433|]
</lang> Further, let's define multiplication and subtraction to check our results... <lang OCaml> let mul m v =
Array.mapi (fun i u -> Array.foldi_left (fun j sum uj -> sum +. uj *. v.(j) ) 0. u ) m
let sub u v = Array.mapi (fun i e -> e -. v.(i)) u </lang> Now 'x' can be plugged into the equation to calculate the residual: <lang OCaml>
- let residual = sub b (mul a x);;
val residual : float array =
[|9.8879238130678e-17; 1.11022302462515654e-16; 2.22044604925031308e-16; 8.88178419700125232e-16; -5.5511151231257827e-16; 4.26741975090294545e-16|]
</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.
Perl
<lang Perl>use Math::Matrix; my $a = Math::Matrix->new([0,1,0],
[0,0,1], [2,0,1]);
my $b = Math::Matrix->new([1],
[2], [4]);
my $x = $a->concat($b)->solve; print $x;</lang>
Math::Matrix
solve()
expects the column vector to be an extra column in the matrix, hence concat()
. Putting not just a column there but a whole identity matrix (making Nx2N) is how its invert()
is implemented. Note that solve()
doesn't notice singular matrices and still gives a return when there is in fact no solution to Ax=B.
Perl 6
Multi-dimensional arrays are not fully implemented in rakudo yet, and that makes Linear Algebra and matrix manipulation quite painful currently. However, it's always possible to proceed as in C with one-dimensional arrays and some index reordering tricks.
<lang perl6>sub mat_elem(@a, $y, $x, $n) is rw { @a[ $y * $n + $x ] } sub swap_row(@a, @b, $r1, $r2, $n) {
return if $r1 == $r2; for ^$n -> $i {
( mat_elem(@a, $r1, $i, $n), mat_elem(@a, $r2, $i, $n) ).=reverse;
} @b[$r1, $r2].=reverse;
}
sub gauss_eliminate(@a, @b, $n) {
sub A($y, $x) is rw { mat_elem(@a, $y, $x, $n) } my ($i, $j, $col, $row, $max_row, $dia); my ($max, $tmp); for ^$n -> $dia {
for $dia ^..^ $n -> $row { swap_row @a, @b, $dia, max(:by({ abs(A($_, $dia)) }), $dia ^..^ $n), $n; $tmp = A($row, $dia) / A($dia, $dia); for $dia ^..^ $n -> $col { A($row, $col) -= $tmp * A($dia, $col); } A($row, $dia) = 0; @b[$row] -= $tmp * @b[$dia]; }
} my @x; for $n - 1, $n - 2 ... 0 -> $row {
$tmp = @b[$row]; for $n - 1, $n - 2 ...^ $row -> $j { $tmp -= @x[$j] * A($row, $j); } @x[$row] = $tmp / A($row, $row);
} return @x;
}
sub MAIN {
my @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 >; my @b = < -0.01 0.61 0.91 0.99 0.60 0.02 >; .say for gauss_eliminate(@a, @b, 6);
}</lang>
- Output:
-0.01 1.60279039450211 -1.61320305990556 1.24549412137144 -0.490989719584658 0.0657606961752
PHP
<lang php>function swap_rows(&$a, &$b, $r1, $r2) {
if ($r1 == $r2) return;
$tmp = $a[$r1]; $a[$r1] = $a[$r2]; $a[$r2] = $tmp;
$tmp = $b[$r1]; $b[$r1] = $b[$r2]; $b[$r2] = $tmp;
}
function gauss_eliminate($A, $b, $N) {
for ($col = 0; $col < $N; $col++) { $j = $col; $max = $A[$j][$j];
for ($i = $col + 1; $i < $N; $i++) { $tmp = abs($A[$i][$col]); if ($tmp > $max) { $j = $i; $max = $tmp; } }
swap_rows($A, $b, $col, $j);
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]; } } $x = array(); 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]; } return $x;
}
function test_gauss() {
$a = array( array(1.00, 0.00, 0.00, 0.00, 0.00, 0.00), array(1.00, 0.63, 0.39, 0.25, 0.16, 0.10), array(1.00, 1.26, 1.58, 1.98, 2.49, 3.13), array(1.00, 1.88, 3.55, 6.70, 12.62, 23.80), array(1.00, 2.51, 6.32, 15.88, 39.90, 100.28), array(1.00, 3.14, 9.87, 31.01, 97.41, 306.02) ); $b = array( -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 );
$x = gauss_eliminate($a, $b, 6);
ksort($x); print_r($x);
}
test_gauss();</lang>
- Output:
Array ( [0] => -0.01 [1] => 1.6027903945021 [2] => -1.6132030599055 [3] => 1.2454941213714 [4] => -0.49098971958463 [5] => 0.065760696175228 )
PL/I
<lang pli>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;</lang>
- Output:
Program to solve n simultaneous equations of the form Ax = b. Please type n: Please type A: Please type the right-hand sides, b: The equations are: 1 2 3 14 2 1 3 13 3 -2 -1 -4 Solutions: X(1)= 1.00000000000000000E+0000 X(2)= 2.00000000000000000E+0000 X(3)= 3.00000000000000000E+0000; Residuals: 0.00000000000000000E+0000 0.00000000000000000E+0000 0.00000000000000000E+0000
PowerShell
<lang PowerShell> function gauss($a,$b) {
$n = $a.count for ($k = 0; $k -lt $n; $k++) { $lmax, $max = $k, [Math]::Abs($a[$k][$k]) for ($l = $k+1; $l -lt $n; $l++) { $tmp = [Math]::Abs($a[$l][$k]) if($max -lt $tmp) { $max, $lmax = $tmp, $l } } if ($k -ne $lmax) { $a[$k], $a[$lmax] = $a[$lmax], $a[$k] $b[$k], $b[$lmax] = $b[$lmax], $b[$k] } $akk = $a[$k][$k] for ($i = $k+1; $i -lt $n; $i++){ $aik = $a[$i][$k] for ($j = $k; $j -lt $n; $j++) { $a[$i][$j] = $a[$i][$j]*$akk - $a[$k][$j]*$aik } $b[$i] = $b[$i]*$akk - $b[$k]*$aik } } for ($i = $n-1; $i -ge 0; $i--) { for ($j = $i+1; $j -le $n; $j++) { $b[$i] -= $b[$j]*$a[$i][$j] } $b[$i] = $b[$i]/$a[$i][$i] } $b
} function show($a) {
if($a) { 0..($a.Count - 1) | foreach{ if($a[$_]){"$($a[$_][0..($a[$_].count -1)])"}else{""} } }
} $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) ) "a =" show $a "" $b = @(-0.01, 0.61, 0.91, 0.99, 0.60, 0.02) "b =" $b "" "x =" gauss $a $b
</lang> Output:
a = 1 0 0 0 0 0 1 0.63 0.39 0.25 0.16 0.1 1 1.26 1.58 1.98 2.49 3.13 1 1.88 3.55 6.7 12.62 23.8 1 2.51 6.32 15.88 39.9 100.28 1 3.14 9.87 31.01 97.41 306.02 b = -0.01 0.61 0.91 0.99 0.6 0.02 x = -0.01 1.60279039450213 -1.6132030599056 1.24549412137148 -0.490989719584674 0.0657606961752342
Python
<lang python># The 'gauss' function takes two matrices, 'a' and 'b', with 'a' square, and it return the determinant of 'a' and a matrix 'x' such that a*x = b.
- If 'b' is the identity, then 'x' is the inverse of 'a'.
import copy from fractions import Fraction
def gauss(a, b):
a = copy.deepcopy(a) b = copy.deepcopy(b) n = len(a) p = len(b[0]) det = 1 for i in range(n - 1): k = i for j in range(i + 1, n): if abs(a[j][i]) > abs(a[k][i]): k = j if k != i: a[i], a[k] = a[k], a[i] b[i], b[k] = b[k], b[i] det = -det for j in range(i + 1, n): t = a[j][i]/a[i][i] for k in range(i + 1, n): a[j][k] -= t*a[i][k] for k in range(p): b[j][k] -= t*b[i][k] for i in range(n - 1, -1, -1): for j in range(i + 1, n): t = a[i][j] for k in range(p): b[i][k] -= t*b[j][k] t = 1/a[i][i] det *= a[i][i] for j in range(p): b[i][j] *= t return det, b
def zeromat(p, q):
return [[0]*q for i in range(p)]
def matmul(a, b):
n, p = len(a), len(a[0]) p1, q = len(b), len(b[0]) if p != p1: raise ValueError("Incompatible dimensions") c = zeromat(n, q) for i in range(n): for j in range(q): c[i][j] = sum(a[i][k]*b[k][j] for k in range(p)) return c
def mapmat(f, a):
return [list(map(f, v)) for v in a]
def ratmat(a):
return mapmat(Fraction, a)
- As an example, compute the determinant and inverse of 3x3 magic square
a = [[2, 9, 4], [7, 5, 3], [6, 1, 8]] b = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] det, c = gauss(a, b)
det -360.0
c [[-0.10277777777777776, 0.18888888888888888, -0.019444444444444438], [0.10555555555555554, 0.02222222222222223, -0.061111111111111116], [0.0638888888888889, -0.14444444444444446, 0.14722222222222223]]
- Check product
matmul(a, c) [[1.0, 0.0, 0.0], [5.551115123125783e-17, 1.0, 0.0], [1.1102230246251565e-16, -2.220446049250313e-16, 1.0]]
- Same with fractions, so the result is exact
det, c = gauss(ratmat(a), ratmat(b))
det Fraction(-360, 1)
c [[Fraction(-37, 360), Fraction(17, 90), Fraction(-7, 360)], [Fraction(19, 180), Fraction(1, 45), Fraction(-11, 180)], [Fraction(23, 360), Fraction(-13, 90), Fraction(53, 360)]]
matmul(a, c) [[Fraction(1, 1), Fraction(0, 1), Fraction(0, 1)], [Fraction(0, 1), Fraction(1, 1), Fraction(0, 1)], [Fraction(0, 1), Fraction(0, 1), Fraction(1, 1)]]</lang>
Racket
<lang racket>
- lang racket
(require math/matrix) (define A
(matrix [[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]]))
(define b (col-matrix [-0.01 0.61 0.91 0.99 0.60 0.02]))
(matrix-solve A b) </lang>
- Output:
<lang racket>
- <array
'#(6 1) #[-0.01 1.602790394502109 -1.613203059905556 1.2454941213714346 -0.4909897195846582 0.06576069617523222]>
</lang>
REXX
version 1
<lang rexx>/* REXX ---------------------------------------------------------------
- 07.08.2014 Walter Pachl translated from PL/I)
- --------------------------------------------------------------------*/
Numeric Digits 20 Parse Arg t n=3 Parse Value '1 2 3 14' With a.1.1 a.1.2 a.1.3 b.1 Parse Value '2 1 3 13' With a.2.1 a.2.2 a.2.3 b.2 Parse Value '3 -2 -1 -4' With a.3.1 a.3.2 a.3.3 b.3 If t=6 Then Do n=6 Parse Value '1.00 0.00 0.00 0.00 0.00 0.00 ' With a.1.1 a.1.2 a.1.3 a.1.4 a.1.5 a.1.6 . Parse Value '1.00 0.63 0.39 0.25 0.16 0.10 ' With a.2.1 a.2.2 a.2.3 a.2.4 a.2.5 a.2.6 . Parse Value '1.00 1.26 1.58 1.98 2.49 3.13 ' With a.3.1 a.3.2 a.3.3 a.3.4 a.3.5 a.3.6 . Parse Value '1.00 1.88 3.55 6.70 12.62 23.80 ' With a.4.1 a.4.2 a.4.3 a.4.4 a.4.5 a.4.6 . Parse Value '1.00 2.51 6.32 15.88 39.90 100.28' With a.5.1 a.5.2 a.5.3 a.5.4 a.5.5 a.5.6 . Parse Value '1.00 3.14 9.87 31.01 97.41 306.02' With a.6.1 a.6.2 a.6.3 a.6.4 a.6.5 a.6.6 . Parse Value '-0.01 0.61 0.91 0.99 0.60 0.02' With b.1 b.2 b.3 b.4 b.5 b.6 . End Do i=1 To n Do j=1 To n sa.i.j=a.i.j End sb.i=b.i End Say 'The equations are:' do i = 1 to n; ol= Do j=1 To n ol=ol format(a.i.j,4,4) End ol=ol' 'format(b.i,4,4) Say ol end
call Gauss_elimination
call Backward_substitution
Say 'Solutions:' Do i=1 To n Say 'x('i')='||x.i End
/* Check solutions: */ Say 'Residuals:' do i = 1 to n res=0 Do j=1 To n res=res+(sa.i.j*x.j) End res=res-sb.i Say 'res('i')='res End
Exit
Gauss_elimination:
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 Return
Backward_substitution:
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 Return</lang>
- Output:
The equations are: 1.0000 2.0000 3.0000 14.0000 2.0000 1.0000 3.0000 13.0000 3.0000 -2.0000 -1.0000 -4.0000 Solutions: x(1)=1 x(2)=2 x(3)=3 Residuals: res(1)=0 res(2)=0 res(3)=0
and with test data from PHP
The equations are: 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0100 1.0000 0.6300 0.3900 0.2500 0.1600 0.1000 0.6100 1.0000 1.2600 1.5800 1.9800 2.4900 3.1300 0.9100 1.0000 1.8800 3.5500 6.7000 12.6200 23.8000 0.9900 1.0000 2.5100 6.3200 15.8800 39.9000 100.2800 0.6000 1.0000 3.1400 9.8700 31.0100 97.4100 306.0200 0.0200 Solutions: x(1)=-0.01 x(2)=1.6027903945021139463 x(3)=-1.6132030599055614262 x(4)=1.2454941213714367527 x(5)=-0.49098971958465761669 x(6)=0.065760696175232005188 Residuals: res(1)=0 res(2)=0.00000000000000000001 res(3)=-0.00000000000000000016 res(4)=0 res(5)=-0.0000000000000000017 res(6)=0.000000000000000001
version 2
Programming note: with the large precision (numeric digits 200), the residuals were insignificant. <lang rexx>/*REXX program solves Ax=b with Gaussian elimination & backwards substitution.*/ parse arg iFID .; if iFID== then iFID='GAUSS_E.DAT' /*¬given? Use default.*/ numeric digits 200 /*use hefty precision.*/
do rec=1 while lines(iFID)\==0 /*read equation sets. */ #=0 do $=1 while lines(iFID)\==0 /*process the equation*/ z=linein(iFID); if z= then leave /*Blank line? e─o─data*/ if $==1 then do; say; say center(' equations ',75,'▒'); say; end say z /*show an equation. */ if left(space(z),1)=='*' then iterate /*ignore any comments.*/ #=#+1; n=words(z)-1 /*assign equation #s. */ do e=1 for n; a.#.e=word(z,e); end /*e*/ /*process the A #s. */ b.#=word(z,n+1) /* " " B #s. */ end /*$*/ if #\==0 then call Gauss_elimination /*compute,show results*/ end /*rec*/
exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ Gauss_elimination: do j=1 for n; jp=j+1
do i=jp to n; _=a.j.j/a.i.j do k=jp to n; a.i.k=a.j.k-_*a.i.k; end /*k*/ b.i=b.j-_*b.i end /*i*/ end /*j*/
x.n=b.n/a.n.n
do j=n-1 to 1 by -1; _=0 do i=j+1 to n; _=_+a.j.i*x.i; end /*i*/ x.j=(b.j-_)/a.j.j end /*j*/ /* [↑] uses backwards substitution. */
numeric digits 8 /*for the display, only use 8 digits.*/ say; say center('solution',75,'═'); say /*a title line for articulated output*/
do o=1 for n; say right('x['o"] = ",38) left(,x.o>=0) x.o/1; end /*o*/
return</lang> input file: GAUSS_E.DAT
* a1 a2 a3 b * ─── ─── ─── ─── 1 2 3 14 2 1 3 13 3 -2 -1 -4 * a1 a2 a3 a4 a5 a6 b * ─────── ─────── ─────── ─────── ─────── ─────── ─────── 1 0 0 0 0 0 -0.01 1 0.63 0.39 0.25 0.16 0.10 0.61 1 1.26 1.58 1.98 2.49 3.13 0.91 1 1.88 3.55 6.70 12.62 23.80 0.99 1 2.51 6.32 15.88 39.90 100.28 0.60 1 3.14 9.87 31.01 97.41 306.02 0.02
output when using the default input file:
▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ equations ▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ * a1 a2 a3 b * ─── ─── ─── ─── 1 2 3 14 2 1 3 13 3 -2 -1 -4 ═════════════════════════════════solution══════════════════════════════════ x[1] = 1 x[2] = 2 x[3] = 3 ▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ equations ▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ * a1 a2 a3 a4 a5 a6 b * ─────── ─────── ─────── ─────── ─────── ─────── ─────── 1 0 0 0 0 0 -0.01 1 0.63 0.39 0.25 0.16 0.10 0.61 1 1.26 1.58 1.98 2.49 3.13 0.91 1 1.88 3.55 6.70 12.62 23.80 0.99 1 2.51 6.32 15.88 39.90 100.28 0.60 1 3.14 9.87 31.01 97.41 306.02 0.02 ═════════════════════════════════solution══════════════════════════════════ x[1] = -0.01 x[2] = 1.6027904 x[3] = -1.6132031 x[4] = 1.2454941 x[5] = -0.49098972 x[6] = 0.065760696
Ruby
<lang ruby> require 'bigdecimal/ludcmp' include LUSolve
BigDecimal::limit(30)
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].map{|i|BigDecimal(i,16)}
b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02].map{|i|BigDecimal(i,16)}
n = 6 zero = BigDecimal("0.0") one = BigDecimal("1.0")
lusolve(a, b, ludecomp(a, n, zero,one), zero).each{|v| puts v.to_s('F')[0..20]}</lang>
- Output:
-0.01 1.6027903945021135753 -1.613203059905560094 1.2454941213714351826 -0.490989719584656871 0.0657606961752318825
Sidef
<lang ruby>var Matrix = require('Math::Matrix');
var a = Matrix.new([0,1,0],
[0,0,1], [2,0,1]);
var b = Matrix.new([1],
[2], [4]);
a.concat(b).solve.print;</lang>
- Output:
1.00000 1.00000 2.00000
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
TI-83 BASIC
The rref() function performs reduced row-echelon form using Gaussian elimination
on a n*(n+1) matrix. The (n+1)th column receives the resulting vector.
The n*n maxtrix is set to 0 and the pivots are set to 1.
The Matr>List() subroutine extracts the (n+1)th column to a list.
The matrix can be more easily entered by the matrix editor.
On TI-83 or TI-84, another way to solve this task is to use the PlySmlt2 internal apps and choose
"simult equ solver" with 6 equations and 6 unknowns.
<lang ti83b>[[ 1.00 0.00 0.00 0.00 0.00 0.00 -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.90 100.28 0.60] [ 1.00 3.14 9.87 31.01 97.41 306.02 0.02]]→[A]
Matr>List(rref([A]),7,L1) L1</lang>
- Output:
{-.01 1.602790395 -1.61320306 1.245494121 -.4909897196 .0657606962}
zkl
<lang zkl>fcn gaussEliminate(a,b){ // modifies a&b --> vector
n:=b.len(); foreach dia in ([0..n-1]){ maxRow:=dia; max:=a[dia][dia]; foreach row in ([dia+1 .. n-1]){ if((tmp:=a[row][dia].abs()) > max){ maxRow=row; max=tmp; } } a.swap(dia,maxRow); b.swap(dia,maxRow); // swap rows foreach row in ([dia+1 .. n-1]){ ar:=a[row]; ad:=a[dia]; tmp:=ar[dia] / ad[dia];
foreach col in ([dia+1 .. n-1]){ ar[col]-=tmp*ad[col]; } ar[dia]=0.0; b[row]-=tmp*b[dia];
} } x:=(0).pump(n,List().write); // -->list filled with garbage foreach row in ([n-1 .. 0,-1]){ tmp:=b[row]; ar:=a[row]; foreach j in ([n-1 .. row+1,-1]){ tmp-=x[j]*ar[j]; } x[row]=tmp/a[row][row]; } x
}</lang> <lang zkl>a:=List( List(1.00, 0.00, 0.00, 0.00, 0.00, 0.00,),
List(1.00, 0.63, 0.39, 0.25, 0.16, 0.10,), List(1.00, 1.26, 1.58, 1.98, 2.49, 3.13,), List(1.00, 1.88, 3.55, 6.70, 12.62, 23.80,), List(1.00, 2.51, 6.32, 15.88, 39.90, 100.28,), List(1.00, 3.14, 9.87, 31.01, 97.41, 306.02) );
b:=List( -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 ); gaussEliminate(a,b).println();</lang>
- Output:
L(-0.01,1.60279,-1.6132,1.24549,-0.49099,0.0657607)