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
D
<lang d>import std.stdio, std.math, std.algorithm, std.range, std.numeric;
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 => cast(double[])(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!((a, b) => a[k] > b[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]
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>
Julia
Using built-in LAPACK-based linear solver (which employs partial-pivoted Gaussian elimination): <lang julia>x = A \ b</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.
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 PL/I>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
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>
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
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