Gaussian elimination: Difference between revisions
Content added Content deleted
(Java) |
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
||
Line 14: | Line 14: | ||
:* the Wikipedia entry: [[wp:Gaussian elimination|<u>Gaussian elimination</u>]] |
:* the Wikipedia entry: [[wp:Gaussian elimination|<u>Gaussian elimination</u>]] |
||
<br><br> |
<br><br> |
||
=={{header| Lobster}}== |
|||
{{trans|Go}} |
|||
<lang Lobster>import std |
|||
// test case from Go version at http://rosettacode.org/wiki/Gaussian_elimination |
|||
// |
|||
let ta = [[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 tb = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02] |
|||
let tx = [-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. |
|||
let ε = 1.0e-14 |
|||
def GaussPartial(a0, b0) -> [float], string: |
|||
// make augmented matrix |
|||
let m = length(b0) |
|||
let a = map(m): [] |
|||
for(a0) ai, i: |
|||
//let ai = a0[i] |
|||
a[i] = map(m+1) j: if j < m: ai[j] else: b0[i] |
|||
// WP algorithm from Gaussian elimination page produces row-eschelon form |
|||
var i = 0 |
|||
var j = 0 |
|||
for(a0) ak, k: |
|||
// Find pivot for column k: |
|||
var iMax = 0 |
|||
var kmax = -1.0 |
|||
i = k |
|||
while i < m: |
|||
let row = a[i] |
|||
// compute scale factor s = max abs in row |
|||
var s = -1.0 |
|||
j = k |
|||
while j < m: |
|||
s = max(s, abs(row[j])) |
|||
j += 1 |
|||
// scale the abs used to pick the pivot |
|||
let kabs = abs(row[k]) / s |
|||
if kabs > kmax: |
|||
iMax = i |
|||
kmax = kabs |
|||
i += 1 |
|||
if a[iMax][k] == 0: |
|||
return [], "singular" |
|||
// swap rows(k, i_max) |
|||
let row = a[k] |
|||
a[k] = a[iMax] |
|||
a[iMax] = row |
|||
// Do for all rows below pivot: |
|||
i = k + 1 |
|||
while i < m: |
|||
// Do for all remaining elements in current row: |
|||
j = k + 1 |
|||
while j <= m: |
|||
a[i][j] -= a[k][j] * (a[i][k] / a[k][k]) |
|||
j += 1 |
|||
// Fill lower triangular matrix with zeros: |
|||
a[i][k] = 0 |
|||
i += 1 |
|||
// end of WP algorithm; now back substitute to get result |
|||
let x = map(m): 0.0 |
|||
i = m - 1 |
|||
while i >= 0: |
|||
x[i] = a[i][m] |
|||
j = i + 1 |
|||
while j < m: |
|||
x[i] -= a[i][j] * x[j] |
|||
j += 1 |
|||
x[i] /= a[i][i] |
|||
i -= 1 |
|||
return x, "" |
|||
def test(): |
|||
let x, err = GaussPartial(ta, tb) |
|||
if err != "": |
|||
print("Error: " + err) |
|||
return |
|||
print(x) |
|||
for(x) xi, i: |
|||
if abs(tx[i]-xi) > ε: |
|||
print("out of tolerance, expected: " + tx[i] + " got: " + xi) |
|||
test()</lang> |
|||
{{out}} |
|||
<pre> |
|||
[-0.01, 1.602790394502, -1.613203059906, 1.245494121371, -0.490989719585, 0.065760696175] |
|||
</pre> |
|||
=={{header|11l}}== |
=={{header|11l}}== |
||
Line 622: | Line 718: | ||
</pre> |
</pre> |
||
=={{header|C sharp|C#}}== |
|||
=={{header|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> |
|||
{{out}} |
|||
<pre> |
|||
(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) |
|||
</pre> |
|||
=={{header|C#}}== |
|||
This modifies A and b in place, which might not be quite desirable. |
This modifies A and b in place, which might not be quite desirable. |
||
<lang csharp> |
<lang csharp> |
||
Line 864: | Line 926: | ||
-0.553874184782179 |
-0.553874184782179 |
||
0.0723048745759396 |
0.0723048745759396 |
||
</pre> |
|||
=={{header|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> |
|||
{{out}} |
|||
<pre> |
|||
(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) |
|||
</pre> |
</pre> |
||
Line 1,127: | Line 1,223: | ||
x[8]=23.00000000000000 |
x[8]=23.00000000000000 |
||
</pre> |
</pre> |
||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
Gaussian Elimination with partial pivoting using augmented matrix |
Gaussian Elimination with partial pivoting using augmented matrix |
||
Line 1,252: | Line 1,349: | ||
</pre> |
</pre> |
||
=={{header|Go}}== |
=={{header|Go}}== |
||
Line 2,274: | Line 2,369: | ||
<pre> |
<pre> |
||
[-0.01, 1.6027903945021138, -1.6132030599055616, 1.2454941213714392, -0.49098971958465953, 0.06576069617523238] |
[-0.01, 1.6027903945021138, -1.6132030599055616, 1.2454941213714392, -0.49098971958465953, 0.06576069617523238] |
||
</pre> |
|||
=={{header| Lobster}}== |
|||
{{trans|Go}} |
|||
<lang Lobster>import std |
|||
// test case from Go version at http://rosettacode.org/wiki/Gaussian_elimination |
|||
// |
|||
let ta = [[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 tb = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02] |
|||
let tx = [-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. |
|||
let ε = 1.0e-14 |
|||
def GaussPartial(a0, b0) -> [float], string: |
|||
// make augmented matrix |
|||
let m = length(b0) |
|||
let a = map(m): [] |
|||
for(a0) ai, i: |
|||
//let ai = a0[i] |
|||
a[i] = map(m+1) j: if j < m: ai[j] else: b0[i] |
|||
// WP algorithm from Gaussian elimination page produces row-eschelon form |
|||
var i = 0 |
|||
var j = 0 |
|||
for(a0) ak, k: |
|||
// Find pivot for column k: |
|||
var iMax = 0 |
|||
var kmax = -1.0 |
|||
i = k |
|||
while i < m: |
|||
let row = a[i] |
|||
// compute scale factor s = max abs in row |
|||
var s = -1.0 |
|||
j = k |
|||
while j < m: |
|||
s = max(s, abs(row[j])) |
|||
j += 1 |
|||
// scale the abs used to pick the pivot |
|||
let kabs = abs(row[k]) / s |
|||
if kabs > kmax: |
|||
iMax = i |
|||
kmax = kabs |
|||
i += 1 |
|||
if a[iMax][k] == 0: |
|||
return [], "singular" |
|||
// swap rows(k, i_max) |
|||
let row = a[k] |
|||
a[k] = a[iMax] |
|||
a[iMax] = row |
|||
// Do for all rows below pivot: |
|||
i = k + 1 |
|||
while i < m: |
|||
// Do for all remaining elements in current row: |
|||
j = k + 1 |
|||
while j <= m: |
|||
a[i][j] -= a[k][j] * (a[i][k] / a[k][k]) |
|||
j += 1 |
|||
// Fill lower triangular matrix with zeros: |
|||
a[i][k] = 0 |
|||
i += 1 |
|||
// end of WP algorithm; now back substitute to get result |
|||
let x = map(m): 0.0 |
|||
i = m - 1 |
|||
while i >= 0: |
|||
x[i] = a[i][m] |
|||
j = i + 1 |
|||
while j < m: |
|||
x[i] -= a[i][j] * x[j] |
|||
j += 1 |
|||
x[i] /= a[i][i] |
|||
i -= 1 |
|||
return x, "" |
|||
def test(): |
|||
let x, err = GaussPartial(ta, tb) |
|||
if err != "": |
|||
print("Error: " + err) |
|||
return |
|||
print(x) |
|||
for(x) xi, i: |
|||
if abs(tx[i]-xi) > ε: |
|||
print("out of tolerance, expected: " + tx[i] + " got: " + xi) |
|||
test()</lang> |
|||
{{out}} |
|||
<pre> |
|||
[-0.01, 1.602790394502, -1.613203059906, 1.245494121371, -0.490989719585, 0.065760696175] |
|||
</pre> |
</pre> |
||
Line 3,259: | Line 3,258: | ||
<code>Math::Matrix</code> <code>solve()</code> expects the column vector to be an extra column in the matrix, hence <code>concat()</code>. Putting not just a column there but a whole identity matrix (making Nx2N) is how its <code>invert()</code> is implemented. Note that <code>solve()</code> doesn't notice singular matrices and still gives a return when there is in fact no solution to Ax=B. |
<code>Math::Matrix</code> <code>solve()</code> expects the column vector to be an extra column in the matrix, hence <code>concat()</code>. Putting not just a column there but a whole identity matrix (making Nx2N) is how its <code>invert()</code> is implemented. Note that <code>solve()</code> doesn't notice singular matrices and still gives a return when there is in fact no solution to Ax=B. |
||
=={{header|Perl 6}}== |
|||
{{works with|Rakudo|2018.03}} |
|||
Gaussian elimination results in a matrix in row echelon form. Gaussian elimination with back-substitution (also known as Gauss-Jordan elimination) results in a matrix in reduced row echelon form. That being the case, we can reuse much of the code from the [[Reduced row echelon form]] task. Perl 6 stores and does calculations on decimal numbers within its limit of precision using Rational numbers by default, meaning the calculations are exact. |
|||
<lang perl6>sub gauss-jordan-solve (@a, @b) { |
|||
@b.kv.map: { @a[$^k].append: $^v }; |
|||
@a.&rref[*]»[*-1]; |
|||
} |
|||
# reduced row echelon form (Gauss-Jordan elimination) |
|||
sub rref (@m) { |
|||
return unless @m; |
|||
my ($lead, $rows, $cols) = 0, +@m, +@m[0]; |
|||
for ^$rows -> $r { |
|||
$lead < $cols or return @m; |
|||
my $i = $r; |
|||
until @m[$i;$lead] { |
|||
++$i == $rows or next; |
|||
$i = $r; |
|||
++$lead == $cols and return @m; |
|||
} |
|||
@m[$i, $r] = @m[$r, $i] if $r != $i; |
|||
my $lv = @m[$r;$lead]; |
|||
@m[$r] »/=» $lv; |
|||
for ^$rows -> $n { |
|||
next if $n == $r; |
|||
@m[$n] »-=» @m[$r] »*» (@m[$n;$lead] // 0); |
|||
} |
|||
++$lead; |
|||
} |
|||
@m |
|||
} |
|||
sub rat-or-int ($num) { |
|||
return $num unless $num ~~ Rat; |
|||
return $num.narrow if $num.narrow.WHAT ~~ Int; |
|||
$num.nude.join: '/'; |
|||
} |
|||
sub say-it ($message, @array, $fmt = " %8s") { |
|||
say "\n$message"; |
|||
$_».&rat-or-int.fmt($fmt).put for @array; |
|||
} |
|||
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-it 'A matrix:', @a, "%6.2f"; |
|||
say-it 'or, A in exact rationals:', @a; |
|||
say-it 'B matrix:', @b, "%6.2f"; |
|||
say-it 'or, B in exact rationals:', @b; |
|||
say-it 'x matrix:', (my @gj = gauss-jordan-solve @a, @b), "%16.12f"; |
|||
say-it 'or, x in exact rationals:', @gj, "%28s"; |
|||
</lang> |
|||
{{out}} |
|||
<pre>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 |
|||
or, A in exact rationals: |
|||
1 0 0 0 0 0 |
|||
1 63/100 39/100 1/4 4/25 1/10 |
|||
1 63/50 79/50 99/50 249/100 313/100 |
|||
1 47/25 71/20 67/10 631/50 119/5 |
|||
1 251/100 158/25 397/25 399/10 2507/25 |
|||
1 157/50 987/100 3101/100 9741/100 15301/50 |
|||
B matrix: |
|||
-0.01 |
|||
0.61 |
|||
0.91 |
|||
0.99 |
|||
0.60 |
|||
0.02 |
|||
or, B in exact rationals: |
|||
-1/100 |
|||
61/100 |
|||
91/100 |
|||
99/100 |
|||
3/5 |
|||
1/50 |
|||
x matrix: |
|||
-0.010000000000 |
|||
1.602790394502 |
|||
-1.613203059906 |
|||
1.245494121371 |
|||
-0.490989719585 |
|||
0.065760696175 |
|||
or, x in exact rationals: |
|||
-1/100 |
|||
655870882787/409205648497 |
|||
-660131804286/409205648497 |
|||
509663229635/409205648497 |
|||
-200915766608/409205648497 |
|||
26909648324/409205648497 |
|||
</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Line 3,964: | Line 3,851: | ||
0.06576069617523222]> |
0.06576069617523222]> |
||
</lang> |
</lang> |
||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
{{works with|Rakudo|2018.03}} |
|||
Gaussian elimination results in a matrix in row echelon form. Gaussian elimination with back-substitution (also known as Gauss-Jordan elimination) results in a matrix in reduced row echelon form. That being the case, we can reuse much of the code from the [[Reduced row echelon form]] task. Perl 6 stores and does calculations on decimal numbers within its limit of precision using Rational numbers by default, meaning the calculations are exact. |
|||
<lang perl6>sub gauss-jordan-solve (@a, @b) { |
|||
@b.kv.map: { @a[$^k].append: $^v }; |
|||
@a.&rref[*]»[*-1]; |
|||
} |
|||
# reduced row echelon form (Gauss-Jordan elimination) |
|||
sub rref (@m) { |
|||
return unless @m; |
|||
my ($lead, $rows, $cols) = 0, +@m, +@m[0]; |
|||
for ^$rows -> $r { |
|||
$lead < $cols or return @m; |
|||
my $i = $r; |
|||
until @m[$i;$lead] { |
|||
++$i == $rows or next; |
|||
$i = $r; |
|||
++$lead == $cols and return @m; |
|||
} |
|||
@m[$i, $r] = @m[$r, $i] if $r != $i; |
|||
my $lv = @m[$r;$lead]; |
|||
@m[$r] »/=» $lv; |
|||
for ^$rows -> $n { |
|||
next if $n == $r; |
|||
@m[$n] »-=» @m[$r] »*» (@m[$n;$lead] // 0); |
|||
} |
|||
++$lead; |
|||
} |
|||
@m |
|||
} |
|||
sub rat-or-int ($num) { |
|||
return $num unless $num ~~ Rat; |
|||
return $num.narrow if $num.narrow.WHAT ~~ Int; |
|||
$num.nude.join: '/'; |
|||
} |
|||
sub say-it ($message, @array, $fmt = " %8s") { |
|||
say "\n$message"; |
|||
$_».&rat-or-int.fmt($fmt).put for @array; |
|||
} |
|||
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-it 'A matrix:', @a, "%6.2f"; |
|||
say-it 'or, A in exact rationals:', @a; |
|||
say-it 'B matrix:', @b, "%6.2f"; |
|||
say-it 'or, B in exact rationals:', @b; |
|||
say-it 'x matrix:', (my @gj = gauss-jordan-solve @a, @b), "%16.12f"; |
|||
say-it 'or, x in exact rationals:', @gj, "%28s"; |
|||
</lang> |
|||
{{out}} |
|||
<pre>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 |
|||
or, A in exact rationals: |
|||
1 0 0 0 0 0 |
|||
1 63/100 39/100 1/4 4/25 1/10 |
|||
1 63/50 79/50 99/50 249/100 313/100 |
|||
1 47/25 71/20 67/10 631/50 119/5 |
|||
1 251/100 158/25 397/25 399/10 2507/25 |
|||
1 157/50 987/100 3101/100 9741/100 15301/50 |
|||
B matrix: |
|||
-0.01 |
|||
0.61 |
|||
0.91 |
|||
0.99 |
|||
0.60 |
|||
0.02 |
|||
or, B in exact rationals: |
|||
-1/100 |
|||
61/100 |
|||
91/100 |
|||
99/100 |
|||
3/5 |
|||
1/50 |
|||
x matrix: |
|||
-0.010000000000 |
|||
1.602790394502 |
|||
-1.613203059906 |
|||
1.245494121371 |
|||
-0.490989719585 |
|||
0.065760696175 |
|||
or, x in exact rationals: |
|||
-1/100 |
|||
655870882787/409205648497 |
|||
-660131804286/409205648497 |
|||
509663229635/409205648497 |
|||
-200915766608/409205648497 |
|||
26909648324/409205648497 |
|||
</pre> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
Line 4,727: | Line 4,727: | ||
End Sub</lang>{{out}} |
End Sub</lang>{{out}} |
||
<pre>(-0.01, 1.60279039450209, -1.61320305990548, 1.24549412137136, -0.490989719584628, 0.065760696175228)</pre> |
<pre>(-0.01, 1.60279039450209, -1.61320305990548, 1.24549412137136, -0.490989719584628, 0.065760696175228)</pre> |
||
=={{header|VBScript}}== |
=={{header|VBScript}}== |
||
<lang vb>' Gaussian elimination - VBScript |
<lang vb>' Gaussian elimination - VBScript |