さまざまな言語で数値計算
Only Do What Only You Can Do
ガウス・ジョルダン法
ガウス・ジョルダン法を使って, 連立一次方程式の解を求める .
考え方はガウスの消去法とほぼ同じ.
例として, 連立一次方程式
を考える.
1. 前進消去
まず(1)式を 9 で割って(1)'とする.
(2)式から(1)'式に 2 を掛けたものを引くと $x$ の項が消える.
(3)式に(1)'式を足すと $x$ の項が消え, (4)式から(1)'式を引くと $x$ の項が消える.
ここまでは, ガウスの消去法と同じ.
ガウスの消去法は,
2番目の式を利用して, 3番目以降の式の $y$ を消す.
3番目の式を利用して, 4番目の式の $z$ を消す...
という方法だったが, ガウス・ジョルダン法は,
2番目の式を利用して, 1番目の式と, 3番目以降の式の $y$ を消す.
3番目の式を利用して, 1番目, 2番目の式と, 4番目の式の $z$ を消す...
というようにして 対角行列を得る.
2. 代入
あとは, 各係数で割れば解が得られる.
VBScript
Option Explicit Private Const N = 4 Private a: a = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) Private b: b = Array(8,17,20,16) 'ピボット選択 pivoting a,b WScript.StdOut.WriteLine "pivoting" WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '前進消去 forward_elimination a,b WScript.StdOut.WriteLine "forward elimination" WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '後退代入 backward_substitution a,b WScript.StdOut.WriteLine "X" disp_vector b '前進消去 Private Sub forward_elimination(ByRef a, ByRef b) Dim pivot, row, col For pivot = 0 To (N - 1) For row = 0 To (N - 1) If row <> pivot Then Dim s: s = a(row)(pivot) / a(pivot)(pivot) For col = pivot To (N - 1) a(row)(col) = a(row)(col) - a(pivot)(col) * s Next b(row) = b(row) - b(pivot) * s End If Next Next End Sub '後退代入 Private Sub backward_substitution(ByVal a, ByRef b) Dim pivot For pivot = 0 To (N -1) b(pivot) = b(pivot) / a(pivot)(pivot) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a, ByRef b) Dim pivot, row, col For pivot = 0 To (N - 1) '各列で 一番値が大きい行を 探す Dim max_row: max_row = pivot Dim max_val: max_val = 0 For row = pivot To (N - 1) If Abs(a(row)(pivot)) > max_val Then '一番値が大きい行 max_val = Abs(a(row)(pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp For col = 0 To (N - 1) tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp Next tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp End If Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub
Z:\>cscript //nologo Z:\1004.vbs pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
JScript
var N = 4 var a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]] var b = [8,17,20,16] // ピボット選択 pivoting(a,b) WScript.StdOut.WriteLine("pivoting") WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 前進消去 forward_elimination(a,b) WScript.StdOut.WriteLine("forward elimination") WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 後退代入 backward_substitution(a,b) WScript.StdOut.WriteLine("X") disp_vector(b) // 前進消去 function forward_elimination(a, b) { for (pivot = 0; pivot < N; pivot++) { for (row = 0; row < N; row++) { if (row == pivot) continue var s = a[row][pivot] / a[pivot][pivot] for (col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s b[row] -= b[pivot] * s } } } // 後退代入 function backward_substitution(a, b) { for (pivot = 0; pivot < N; pivot++) b[pivot] /= a[pivot][pivot] } // ピボット選択 function pivoting(a, b) { for (pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val = 0 for (row = pivot; row < N; row++) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp for (col = 0; col < N; col++) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 状態を確認 function disp_progress(a, b) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + a[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine((" " + b[i].toFixed(10)).slice(-14)) } WScript.StdOut.WriteLine() } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 2次元配列を表示 function disp_matrix(matrix) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + matrix[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }
Z:\>cscript //nologo Z:\1004.js pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
PowerShell
set-variable -name N -value 3 -option constant # 前進消去 function forward_elimination($a, $b) { foreach ($pivot in 0..$N) { foreach ($row in 0..$N) { if ($row -eq $pivot) { continue } $s = $a[$row][$pivot] / $a[$pivot][$pivot] foreach ($col in $pivot..$N) { $a[$row][$col] -= $a[$pivot][$col] * $s } $b[$row] -= $b[$pivot] * $s } } } # 後退代入 function backward_substitution($a, $b) { foreach ($pivot in 0..$N) { $b[$pivot] /= $a[$pivot][$pivot] } } # ピボット選択 function pivoting($a, $b) { foreach ($pivot in 0..$N) { # 各列で 一番値が大きい行を 探す $max_row = $pivot $max_val = 0 foreach ($row in $pivot..$N) { if ([Math]::Abs($a[$row][$pivot]) -gt $max_val) { # 一番値が大きい行 $max_val = [Math]::Abs($a[$row][$pivot]) $max_row = $row } } # 一番値が大きい行と入れ替え if ($max_row -ne $pivot) { foreach ($col in 0..$N) { $tmp = $a[$max_row][$col] $a[$max_row][$col] = $a[$pivot][$col] $a[$pivot][$col] = $tmp } $tmp = $b[$max_row] $b[$max_row] = $b[$pivot] $b[$pivot] = $tmp } } } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } $a = (-1.0,-2.0,7.0,-2.0), (1.0,-1.0,-2.0,6.0), (9.0,2.0,1.0,1.0), (2.0,8.0,-2.0,1.0) $b = 8.0, 17.0, 20.0, 16.0 # ピボット選択 pivoting $a $b Write-Host "pivoting" Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 前進消去 forward_elimination $a $b Write-Host "forward elimination" Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 後退代入 backward_substitution $a $b Write-Host "X" disp_vector $b
Z:\>powershell -file Z:\1004.ps1 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Perl
use constant N => 3; my @a = ([-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]); my @b = (8,17,20,16); # ピボット選択 pivoting(\@a,\@b); print "pivoting\n"; print "A\n"; disp_matrix(\@a); print "B\n"; disp_vector(\@b); print "\n"; # 前進消去 forward_elimination(\@a,\@b); print "forward elimination\n"; print "A\n"; disp_matrix(\@a); print "B\n"; disp_vector(\@b); print "\n"; # 後退代入 backward_substitution(\@a,\@b); print "X\n"; disp_vector(\@b); # 前進消去 sub forward_elimination { my ($a, $b) = @_; for $pivot (0..N) { for $row (0..N) { next if ($row == $pivot); my $s = $$a[$row][$pivot] / $$a[$pivot][$pivot]; for $col ($pivot..N) { $$a[$row][$col] -= $$a[$pivot][$col] * $s; } $$b[$row] -= $$b[$pivot] * $s; } } } # 後退代入 sub backward_substitution { my ($a, $b) = @_; for $pivot (0..N) { $$b[$pivot] /= $$a[$pivot][$pivot]; } } # ピボット選択 sub pivoting { my ($a, $b) = @_; for $pivot (0..N) { # 各列で 一番値が大きい行を 探す my $max_row = $pivot; my $max_val = 0; for $row ($pivot..N) { if (abs($$a[$row][$pivot]) > $max_val) { # 一番値が大きい行 $max_val = abs($$a[$row][$pivot]); $max_row = $row; } } # 一番値が大きい行と入れ替え if ($max_row != $pivot) { my $tmp; for $col (0..N) { $tmp = $$a[$max_row][$col]; $$a[$max_row][$col] = $$a[$pivot][$col]; $$a[$pivot][$col] = $tmp; } $tmp = $$b[$max_row]; $$b[$max_row] = $$b[$pivot]; $$b[$pivot] = $tmp; } } } # 2次元配列を表示 sub disp_matrix { my ($matrix) = @_; foreach $row (@$matrix) { foreach $col (@$row) { printf("%14.10f\t", $col); } print "\n"; } } # 1次元配列を表示 sub disp_vector { my ($row) = @_; foreach $col (@$row) { printf("%14.10f\t", $col); } print "\n"; }
Z:\>perl Z:\1004.pl pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
PHP
<?php define("N", 3); $a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]]; $b = [8,17,20,16]; # ピボット選択 pivoting($a,$b); print "pivoting\n"; print "A\n"; disp_matrix($a); print "B\n"; disp_vector($b); print "\n"; # 前進消去 forward_elimination($a,$b); print "forward elimination\n"; print "A\n"; disp_matrix($a); print "B\n"; disp_vector($b); print "\n"; # 後退代入 backward_substitution($a,$b); print "X\n"; disp_vector($b); # 前進消去 function forward_elimination(&$a, &$b) { foreach (range(0, N) as $pivot) { foreach (range(0, N) as $row) { if ($row == $pivot) continue; $s = $a[$row][$pivot] / $a[$pivot][$pivot]; foreach (range($pivot, N) as $col) $a[$row][$col] -= $a[$pivot][$col] * $s; $b[$row] -= $b[$pivot] * $s; } } } # 後退代入 function backward_substitution(&$a, &$b) { foreach (range(0, N) as $pivot) $b[$pivot] /= $a[$pivot][$pivot]; } # ピボット選択 function pivoting(&$a, &$b) { foreach (range(0, N) as $pivot) { # 各列で 一番値が大きい行を 探す $max_row = $pivot; $max_val = 0; foreach (range($pivot, N) as $row) { if (abs($a[$row][$pivot]) > $max_val) { # 一番値が大きい行 $max_val = abs($a[$row][$pivot]); $max_row = $row; } } # 一番値が大きい行と入れ替え if ($max_row != $pivot) { foreach (range(0, N) as $col) { $tmp = $a[$max_row][$col]; $a[$max_row][$col] = $a[$pivot][$col]; $a[$pivot][$col] = $tmp; } $tmp = $b[$max_row]; $b[$max_row] = $b[$pivot]; $b[$pivot] = $tmp; } } } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($matrix as $row) { foreach ($row as $col) printf("%14.10f\t", $col); print "\n"; } } # 1次元配列を表示 function disp_vector($row) { foreach ($row as $col) printf("%14.10f\t", $col); print "\n"; } ?>
Z:\>php Z:\1004.php pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Python
# coding: Shift_JIS N = 4 # 前進消去 def forward_elimination(a, b): for pivot in range (0, N, 1): for row in range (0, N, 1): if (row == pivot): continue s = a[row][pivot] / a[pivot][pivot] for col in range(pivot, N, 1): a[row][col] -= a[pivot][col] * s b[row] -= b[pivot] * s # 後退代入 def backward_substitution(a, b): for pivot in range(0, N, 1): b[pivot] /= a[pivot][pivot] # ピボット選択 def pivoting(a, b): for pivot in range(0, N, 1): # 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0.0 for row in range(pivot, N, 1): if (abs(a[row][pivot]) > max_val): # 一番値が大きい行 max_val = abs(a[row][pivot]) max_row = row # 一番値が大きい行と入れ替え if (max_row != pivot): tmp = 0.0 for col in range(0, N, 1): tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp # 2次元配列を表示 def disp_matrix(matrix): for row in matrix: for col in row: print "%14.10f\t" % col, print "" # 1次元配列を表示 def disp_vector(row): for col in row: print "%14.10f\t" % col, print "" a = [[-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0], [ 9.0, 2.0, 1.0, 1.0], [2.0, 8.0, -2.0, 1.0], ] b = [ 8.0, 17.0, 20.0, 16.0] c = [ 0.0, 0.0, 0.0, 0.0] # ピボット選択 pivoting(a,b) print "pivoting" print "A" disp_matrix(a) print "B" disp_vector(b) print "" # 前進消去 forward_elimination(a,b) print "forward_elimination" print "A" disp_matrix(a) print "B" disp_vector(b) print "" # 後退代入 backward_substitution(a,b) print "X" disp_vector(b)
Z:\>python Z:\1004.py pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward_elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Ruby
# coding: Shift_JIS N = 3 # 2次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 前進消去 def forward_elimination(a, b) 0.upto(N) do |pivot| 0.upto(N) do |row| next if (row == pivot) s = a[row][pivot] / a[pivot][pivot] pivot.upto(N) do |col| a[row][col] -= a[pivot][col] * s end b[row] -= b[pivot] * s end end end # 後退代入 def backward_substitution(a, b) 0.upto(N) do |pivot| b[pivot] /= a[pivot][pivot] end end # ピボット選択 def pivoting(a, b) (0..N).each do |pivot| # 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 (pivot..N).each do |row| if ((a[row][pivot]).abs > max_val) # 一番値が大きい行 max_val = (a[row][pivot]).abs max_row = row end end # 一番値が大きい行と入れ替え if (max_row != pivot) tmp = 0 (0..N).each do |col| tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp end tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp end end end # 状態を確認 def disp_progress(a, b) (0..N).each do |i| (0..N).each do |j| printf("%14.10f\t", a[i][j]) end printf("%14.10f\n", b[i]) end puts "" end a = [[-1.0,-2.0,7.0,-2.0],[1.0,-1.0,-2.0,6.0],[9.0,2.0,1.0,1.0],[2.0,8.0,-2.0,1.0]] b = [8.0,17.0,20.0,16.0] # ピボット選択 pivoting(a,b) puts "pivoting" puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 前進消去 forward_elimination(a,b) puts "forward elimination" puts "A" disp_matrix(a) puts "B" disp_vector(b) puts "" # 後退代入 backward_substitution(a,b) puts "X" disp_vector(b)
Z:\>ruby Z:\1004.rb pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Groovy
N = 3 def a = [[-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0], [ 9.0, 2.0, 1.0, 1.0], [2.0, 8.0, -2.0, 1.0]] as double[][] def b = [ 8.0, 17.0, 20.0, 16.0] as double[] def c = [ 0.0, 0.0, 0.0, 0.0] as double[] // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("forward elimination") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 後退代入 backward_substitution(a,b) println("X") disp_vector(b) // 前進消去 def forward_elimination(a, b) { for (pivot in 0..N) { for (row in 0..N) { if (row == pivot) continue s = a[row][pivot] / a[pivot][pivot] for (col in pivot..N) a[row][col] -= a[pivot][col] * s b[row] -= b[pivot] * s } } } // 後退代入 def backward_substitution(a, b) { for (pivot in 0..N) b[pivot] /= a[pivot][pivot] } // ピボット選択 def pivoting(a, b) { for (pivot in 0..N) { // 各列で 一番値が大きい行を 探す max_row = pivot max_val = 0 for (row in pivot..N) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { tmp = 0 for (col in 0..N) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 1次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf ("%14.10f\t" , col) } println() } } // 2次元配列を表示 def disp_vector(row) { for (col in row) { printf ("%14.10f\t" , col) } println() }
Z:\>groovy Groovy1004.groovy pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Pascal
program Pas1004(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 2次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var i,j:Integer; begin for i := Low(matrix) to High(matrix) do begin for j := Low(matrix) to High(matrix) do write(format('%14.10f'#9, [matrix[i,j]])); writeln(); end; end; // 前進消去 procedure forward_elimination(var a:TwoDimArray; var b:array of Double); var s:Double; pivot, row, col:Integer; begin for pivot := Low(a) to High(a) do begin for row := Low(a) to High(a) do begin if row = pivot then continue; s := a[row,pivot] / a[pivot,pivot]; for col := pivot to High(a) do a[row,col] := a[row,col] - a[pivot,col] * s; b[row] := b[row] - b[pivot] * s; end; end; end; // 後退代入 procedure backward_substitution(a:TwoDimArray; var b:array of Double); var pivot:Integer; begin for pivot := Low(a) to High(a) do b[pivot] := b[pivot] / a[pivot,pivot]; end; // ピボット選択 procedure pivoting(var a:TwoDimArray; var b:array of Double); var max_val, tmp:Double; pivot, row, col, max_row:Integer; begin for pivot := Low(a) to High(a) do begin // 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0; for row := pivot to High(a) do begin if Abs(a[row,pivot]) > max_val then begin // 一番値が大きい行 max_val := Abs(a[row,pivot]); max_row := row; end; end; // 一番値が大きい行と入れ替え if max_row <> pivot then begin tmp := 0; for col := Low(a) to High(a) do begin tmp := a[max_row,col]; a[max_row,col] := a[pivot,col]; a[pivot,col] := tmp; end; tmp := b[max_row]; b[max_row] := b[pivot]; b[pivot] := tmp; end; end; end; var a :TwoDimArray = (( 9.0, 2.0, 1.0, 1.0), (2.0, 8.0, -2.0, 1.0), (-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0)); b :array [0..N] of Double = (20.0, 16.0, 8.0, 17.0); begin // ピボット選択 pivoting(a,b); writeln('pivoting'); writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 前進消去 forward_elimination(a,b); writeln('forward elimination'); writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 後退代入 backward_substitution(a,b); writeln('X'); disp_vector(b); end.
Z:\>fpc -v0 -l- Pas1004.pp Z:\>Pas1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Ada
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1004 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_TwoDimArray := ((-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0), ( 9.0, 2.0, 1.0, 1.0), (2.0, 8.0, -2.0, 1.0)); b:Long_Float_Array := (8.0, 17.0, 20.0, 16.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 2次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for i in matrix'Range loop for j in matrix'Range loop Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- 前進消去 procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is s:Long_Float; begin for pivot in a'Range loop for row in a'Range loop if row /= pivot then s := a(row,pivot) / a(pivot,pivot); for col in pivot..a'Last loop a(row,col) := a(row,col) - a(pivot,col) * s; end loop; b(row) := b(row) - b(pivot) * s; end if; end loop; end loop; end forward_elimination; -- 後退代入 procedure backward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array) is begin for pivot in a'Range loop b(pivot) := b(pivot) / a(pivot,pivot); end loop; end backward_substitution; -- ピボット選択 procedure pivoting(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is max_val, tmp:Long_Float; max_row:Integer; begin for pivot in a'Range loop -- 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0.0; for row in pivot..a'Last loop if Abs(a(row,pivot)) > max_val then -- 一番値が大きい行 max_val := Abs(a(row,pivot)); max_row := row; end if; end loop; -- 一番値が大きい行と入れ替え if max_row /= pivot then tmp := 0.0; for col in a'Range loop tmp := a(max_row,col); a(max_row,col) := a(pivot,col); a(pivot,col) := tmp; end loop; tmp := b(max_row); b(max_row) := b(pivot); b(pivot) := tmp; end if; end loop; end pivoting; begin -- ピボット選択 pivoting(a,b); Put_Line("pivoting"); Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 前進消去 forward_elimination(a,b); Put_Line("forward elimination"); Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 後退代入 backward_substitution(a,b); Put_Line("X"); disp_vector(b); end Ada1004;
xxxxxx@yyyyyy /Z $ gnatmake Ada1004.adb xxxxxx@yyyyyy /Z $ Ada1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 -0.0000000000 0.0000000000 0.0000000000 7.5555555556 -0.0000000000 -0.0000000000 -0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 -0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
VB.NET
Option Explicit Module VB1004 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} Dim b() As Double = {8,17,20,16} 'ピボット選択 pivoting(a,b) Console.WriteLine("pivoting") Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '前進消去 forward_elimination(a,b) Console.WriteLine("forward elimination") Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '後退代入 backward_substitution(a,b) Console.WriteLine("X") disp_vector(b) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '2次元配列の表示 Private Sub disp_matrix(ByVal matrix(,) As Double) Dim i As Integer = 0 For Each col As Double In matrix Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) i += 1 If i > N Then i = 0 Console.WriteLine() End If Next End Sub '前進消去 Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N For row As Integer = 0 To N If row = pivot Then Continue For Dim s As Double = a(row,pivot) / a(pivot,pivot) For col As Integer = pivot To N a(row,col) -= a(pivot,col) * s Next b(row) -= b(pivot) * s Next Next End Sub '後退代入 Private Sub backward_substitution(ByVal a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N b(pivot) /= a(pivot,pivot) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N '各列で 一番値が大きい行を 探す Dim max_row As Integer = pivot Dim max_val As Double = 0 For row As Integer = pivot To N If Math.Abs(a(row,pivot)) > max_val Then '一番値が大きい行 max_val = Math.Abs(a(row,pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp As Double For col As Integer = 0 To N tmp = a(max_row,col) a(max_row,col) = a(pivot,col) a(pivot,col) = tmp Next tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp End If Next End Sub End Module
Z:\>vbc -nologo VB1004.vb Z:\>VB1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
C#
using System; public class CS1004 { private const int N = 4; // ガウス・ジョルダン法 public static void Main() { double[,] a = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double[] b = {8,17,20,16}; // ピボット選択 pivoting(a,b); Console.WriteLine("pivoting"); Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 前進消去 forward_elimination(a,b); Console.WriteLine("forward elimination"); Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 後退代入 backward_substitution(a,b); Console.WriteLine("X"); disp_vector(b); } // 前進消去 private static void forward_elimination(double[,] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) { for (int row = 0; row < N; row++) { if (row == pivot) continue; double s = a[row,pivot] / a[pivot,pivot]; for (int col = pivot; col < N; col++) a[row,col] -= a[pivot,col] * s; b[row] -= b[pivot] * s; } } } // 後退代入 private static void backward_substitution(double[,] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) b[pivot] /= a[pivot,pivot]; } // ピボット選択 private static void pivoting(double[,] a, double[] b) { for(int pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; for (int row = pivot; row < N; row++) { if (Math.Abs(a[row,pivot]) > max_val) { // 一番値が大きい行 max_val = Math.Abs(a[row,pivot]); max_row = row; } } // 一番値が大きい行と入れ替え if (max_row != pivot) { double tmp; for (int col = 0; col < N; col++) { tmp = a[max_row,col]; a[max_row,col] = a[pivot,col]; a[pivot,col] = tmp; } tmp = b[max_row]; b[max_row] = b[pivot]; b[pivot] = tmp; } } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 2次元配列を表示 private static void disp_matrix(double[,] matrix) { int i = 0; foreach (double col in matrix) { Console.Write(string.Format("{0,14:F10}\t", col)); if (++i >= N) { i = 0; Console.WriteLine(); } } } }
Z:\>csc -nologo CS1004.cs Z:\>CS1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Java
import java.lang.*; public class Java1004 { private static final int N = 4; // ガウス・ジョルダン法 public static void main(String []args) { double[][] a = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double[] b = {8,17,20,16}; // ピボット選択 pivoting(a,b); System.out.println("pivoting"); System.out.println("A"); disp_matrix(a); System.out.println("B"); disp_vector(b); System.out.println(); // 前進消去 forward_elimination(a,b); System.out.println("forward elimination"); System.out.println("A"); disp_matrix(a); System.out.println("B"); disp_vector(b); System.out.println(); // 後退代入 backward_substitution(a,b); System.out.println("X"); disp_vector(b); } // 前進消去 private static void forward_elimination(double[][] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) { for (int row = 0; row < N; row++) { if (row == pivot) continue; double s = a[row][pivot] / a[pivot][pivot]; for (int col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s; b[row] -= b[pivot] * s; } } } // 後退代入 private static void backward_substitution(double[][] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) b[pivot] /= a[pivot][pivot]; } // ピボット選択 private static void pivoting(double[][] a, double[] b) { for(int pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; for (int row = pivot; row < N; row++) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]); max_row = row; } } // 一番値が大きい行と入れ替え if (max_row != pivot) { double tmp; for (int col = 0; col < N; col++) { tmp = a[max_row][col]; a[max_row][col] = a[pivot][col]; a[pivot][col] = tmp; } tmp = b[max_row]; b[max_row] = b[pivot]; b[pivot] = tmp; } } } // 1次元配列を表示 private static void disp_vector(double[] row) { for (double col: row) System.out.print(String.format("%14.10f\t", col)); System.out.println(); } // 2次元配列を表示 private static void disp_matrix(double[][] matrix) { for (double[] row: matrix) { for (double col: row) System.out.print(String.format("%14.10f\t", col)); System.out.println(); } } }
Z:\>javac Java1004.java Z:\>java Java1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
C++
#include <iostream> #include <iomanip> #include <math.h> using namespace std; const int N = 4; // ピボット選択 void pivoting(double a[N][N], double b[N]); // 前進消去 void forward_elimination(double a[N][N], double b[N]); // 後退代入 void backward_substitution(double a[N][N], double b[N]); // 1次元配列を表示 void disp_vector(double row[N]); // 2次元配列を表示 void disp_matrix(double matrix[N][N]); // ガウス・ジョルダン法 int main() { double a[N][N] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double b[N] = {8,17,20,16}; // ピボット選択 pivoting(a,b); cout << "pivoting" << endl; cout << "A" << endl; disp_matrix(a); cout << "B" << endl; disp_vector(b); cout << endl; // 前進消去 forward_elimination(a,b); cout << "forward elimination" << endl; cout << "A" << endl; disp_matrix(a); cout << "B" << endl; disp_vector(b); cout << endl; // 後退代入 backward_substitution(a,b); cout << "X" << endl; disp_vector(b); return 0; } // 前進消去 void forward_elimination(double a[N][N], double b[N]) { for (int pivot = 0; pivot < N; pivot++) { for (int row = 0; row < N; row++) { if (row == pivot) continue; double s = a[row][pivot] / a[pivot][pivot]; for (int col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s; b[row] -= b[pivot] * s; } } } // 後退代入 void backward_substitution(double a[N][N], double b[N]) { for (int pivot = 0; pivot < N; pivot++) b[pivot] /= a[pivot][pivot]; } // ピボット選択 void pivoting(double a[N][N], double b[N]) { for(int pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; for (int row = pivot; row < N; row++) { if (fabs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = fabs(a[row][pivot]); max_row = row; } } // 一番値が大きい行と入れ替え if (max_row != pivot) { double tmp; for (int col = 0; col < N; col++) { tmp = a[max_row][col]; a[max_row][col] = a[pivot][col]; a[pivot][col] = tmp; } tmp = b[max_row]; b[max_row] = b[pivot]; b[pivot] = tmp; } } } // 1次元配列を表示 void disp_vector(double row[N]) { for (int i = 0; i < N; i++) cout << setw(14) << fixed << setprecision(10) << row[i] << "\t"; cout << endl; } // 2次元配列を表示 void disp_matrix(double matrix[N][N]) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) cout << setw(14) << fixed << setprecision(10) << matrix[i][j] << "\t"; cout << endl; } }
Z:\>bcc32 -q CP1004.cpp cp1004.cpp: Z:\>CP1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 -0.0000000000 0.0000000000 0.0000000000 7.5555555556 -0.0000000000 -0.0000000000 -0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 -0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Objective-C
#import <Foundation/Foundation.h> #import <math.h> const int N = 4; // ピボット選択 void pivoting(double a[N][N], double b[N]); // 前進消去 void forward_elimination(double a[N][N], double b[N]); // 後退代入 void backward_substitution(double a[N][N], double b[N]); // 1次元配列を表示 void disp_vector(double row[N]); // 2次元配列を表示 void disp_matrix(double matrix[N][N]); // ガウス・ジョルダン法 int main() { double a[4][4] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double b[4] = {8,17,20,16}; // ピボット選択 pivoting(a,b); printf("pivoting\n"); printf("A\n"); disp_matrix(a); printf("B\n"); disp_vector(b); printf("\n"); // 前進消去 forward_elimination(a,b); printf("forward elimination\n"); printf("A\n"); disp_matrix(a); printf("B\n"); disp_vector(b); printf("\n"); // 後退代入 backward_substitution(a,b); printf("X\n"); disp_vector(b); return 0; } // 前進消去 void forward_elimination(double a[N][N], double b[N]) { int pivot, row, col; for (pivot = 0; pivot < N; pivot++) { for (row = 0; row < N; row++) { if (row == pivot) continue; double s = a[row][pivot] / a[pivot][pivot]; for (col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s; b[row] -= b[pivot] * s; } } } // 後退代入 void backward_substitution(double a[N][N], double b[N]) { int pivot; for (pivot = 0; pivot < N; pivot++) b[pivot] /= a[pivot][pivot]; } // ピボット選択 void pivoting(double a[N][N], double b[N]) { int pivot, row, col; for(pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; for (row = pivot; row < N; row++) { if (fabs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = fabs(a[row][pivot]); max_row = row; } } // 一番値が大きい行と入れ替え if (max_row != pivot) { double tmp; for (col = 0; col < N; col++) { tmp = a[max_row][col]; a[max_row][col] = a[pivot][col]; a[pivot][col] = tmp; } tmp = b[max_row]; b[max_row] = b[pivot]; b[pivot] = tmp; } } } // 1次元配列を表示 void disp_vector(double row[N]) { int i; for (i = 0; i < N; i++) { printf("%14.10f\t", row[i]); } printf("\n"); } // 2次元配列を表示 void disp_matrix(double matrix[N][N]) { int i, j; for (i = 0; i < N; i++) { for (j = 0; j < N; j++) printf("%14.10f\t", matrix[i][j]); printf("\n"); } }
xxxxxx@yyyyyy /Z $ gcc -o OC1004 OC1004.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 -0.0000000000 0.0000000000 0.0000000000 7.5555555556 -0.0000000000 -0.0000000000 -0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 -0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
D
import std.stdio; import std.math; const int N = 4; // ガウス・ジョルダン法 void main(string[] args) { double[N][N] a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]]; double[N] b = [8,17,20,16]; // ピボット選択 pivoting(a,b); writefln("pivoting"); writefln("A"); disp_matrix(a); writefln("B"); disp_vector(b); writefln(""); // 前進消去 forward_elimination(a,b); writefln("forward elimination"); writefln("A"); disp_matrix(a); writefln("B"); disp_vector(b); writefln(""); // 後退代入 backward_substitution(a,b); writefln("X"); disp_vector(b); } // 前進消去 void forward_elimination(ref double[N][N] a, ref double[N] b) { foreach (pivot; 0..N) { foreach (row; 0..N) { if (row == pivot) continue; double s = a[row][pivot] / a[pivot][pivot]; foreach (col; pivot..N) a[row][col] -= a[pivot][col] * s; b[row] -= b[pivot] * s; } } } // 後退代入 void backward_substitution(double[N][N] a, ref double[N] b) { foreach (pivot; 0..N) b[pivot] /= a[pivot][pivot]; } // ピボット選択 void pivoting(ref double[N][N] a, ref double[N] b) { foreach (pivot; 0..N) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; foreach (row; pivot..N) { if (fabs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = fabs(a[row][pivot]); max_row = row; } } // 一番値が大きい行と入れ替え if (max_row != pivot) { double tmp; foreach (col; 0..N) { tmp = a[max_row][col]; a[max_row][col] = a[pivot][col]; a[pivot][col] = tmp; } tmp = b[max_row]; b[max_row] = b[pivot]; b[pivot] = tmp; } } } // 1次元配列を表示 void disp_vector(double[] row) { foreach(col; row) writef("%14.10f\t", col); writefln(""); } // 2次元配列を表示 void disp_matrix(double[N][N] matrix) { foreach(row; matrix) { foreach(col; row) writef("%14.10f\t", col); writefln(""); } }
Z:\>dmd D1004.d Z:\>D1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Go
package main import "fmt" import "math" const N = 4 // ガウス・ジョルダン法 func main() { var a [N][N]float64 = [N][N]float64{{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} var b []float64 = []float64{8,17,20,16} // ピボット選択 pivoting(&a,b) fmt.Println("pivoting") fmt.Println("A") disp_matrix(a); fmt.Println("B") disp_vector(b); fmt.Println() // 前進消去 forward_elimination(&a,b) fmt.Println("forward elimination") fmt.Println("A") disp_matrix(a); fmt.Println("B") disp_vector(b); fmt.Println() // 後退代入 backward_substitution(&a,b) fmt.Println("X") disp_vector(b); } // 前進消去 func forward_elimination(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { for row := 0; row < N; row++ { if (row == pivot) { continue } var s = a[row][pivot] / a[pivot][pivot] for col := pivot; col < N; col++ { a[row][col] -= a[pivot][col] * s } b[row] -= b[pivot] * s } } } // 後退代入 func backward_substitution(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { b[pivot] /= a[pivot][pivot] } } // ピボット選択 func pivoting(a *[N][N]float64, b []float64) { for pivot := 0; pivot < N; pivot++ { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val float64 = 0.0 for row := pivot; row < N; row++ { if math.Fabs(a[row][pivot]) > max_val { // 一番値が大きい行 max_val = math.Fabs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp float64 = 0.0 for col := 0; col < N; col++ { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 2次元配列を表示 func disp_matrix(matrix[N][N]float64) { for _, row := range matrix { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } }
Z:\>8g GO1004.go Z:\>8l -o GO1004.exe GO1004.8 Z:\>GO1004 pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Scala
object Scala1004 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) var b:Array[Double] = Array(8,17,20,16) // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("forward elimination") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 後退代入 backward_substitution(a,b) println("X") disp_vector(b) } // 前進消去 def forward_elimination(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { for (row <- 0 to N) { if (row != pivot) { var s:Double = a(row)(pivot) / a(pivot)(pivot) for (col <- pivot to N) a(row)(col) -= a(pivot)(col) * s b(row) -= b(pivot) * s } } } } // 後退代入 def backward_substitution(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) b(pivot) /= a(pivot)(pivot) } // ピボット選択 def pivoting(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { // 各列で 一番値が大きい行を 探す var max_row:Integer = pivot var max_val:Double = 0 for (row <- pivot to N) { if (Math.abs(a(row)(pivot)) > max_val) { // 一番値が大きい行 max_val = Math.abs(a(row)(pivot)) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp:Double = 0 for (col <- 0 to N) { tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp } tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp } } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 2次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) = { matrix.foreach { row => row.foreach { col => print("%14.10f\t".format(col)) } println() } } }
Z:\>scala Scala1004.scala pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
F#
module Fs1004 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 2次元配列を表示 let disp_matrix (matrix:float[][]) = matrix |> Array.iter (fun row -> row |> Array.iter (fun col -> printf "%14.10f" col) printfn "" ) // 前進消去 let forward_elimination (a:float[][]) (b:float[]) = for pivot in [0..N] do for row in [0..N] do if row <> pivot then let s:float = a.[row].[pivot] / a.[pivot].[pivot] for col in [pivot..N] do a.[row].[col] <- a.[row].[col] - a.[pivot].[col] * s b.[row] <- b.[row] - b.[pivot] * s // 後退代入 let backward_substitution (a:float[][]) (b:float[]) = for pivot in [0..N] do b.[pivot] <- b.[pivot] / a.[pivot].[pivot] // ピボット選択 let pivoting (a:float[][]) (b:float[]) = for pivot in [0..N] do // 各列で 一番値が大きい行を 探す let mutable max_row:int = pivot let mutable max_val:float = 0.0 for row in [pivot..N] do if (Math.Abs(a.[row].[pivot]) > max_val) then // 一番値が大きい行 max_val <- Math.Abs(a.[row].[pivot]) max_row <- row // 一番値が大きい行と入れ替え if (max_row <> pivot) then let mutable tmp:float = 0.0 for col in [0..N] do tmp <- a.[max_row].[col] a.[max_row].[col] <- a.[pivot].[col] a.[pivot].[col] <- tmp tmp <- b.[max_row] b.[max_row] <- b.[pivot] b.[pivot] <- tmp let mutable a:float[][] = [| [|-1.0;-2.0;7.0;-2.0 |]; [|1.0;-1.0;-2.0;6.0 |]; [| 9.0;2.0;1.0;1.0 |]; [|2.0;8.0;-2.0;1.0 |] |] let mutable b:float[] = [| 8.0;17.0;20.0;16.0 |] // ピボット選択 pivoting a b printfn "pivoting" printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 前進消去 forward_elimination a b printfn "forward elimination" printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 後退代入 backward_substitution a b printfn "X" disp_vector b exit 0
Z:\>fsi --nologo --quiet Fs1004.fs pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Clojure
(def N 4) (def a [[-1.0 -2.0 7.0 -2.0] [1.0 -1.0 -2.0 6.0] [9.0 2.0 1.0 1.0] [2.0 8.0 -2.0 1.0]]) (def b [8.0 17.0 20.0 16.0]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (defn disp_matrix [matrix] (doseq [row matrix] (doseq [col row] (print (format "%14.10f\t" col))) (println))) ;各列で 一番値が大きい行を 探す (defn get_max_row [row col a max_row max_val] ;一番値が大きい行 (def ws1 (if (< max_val (. Math abs (nth (nth a row) col))) (vector row (. Math abs (nth (nth a row) col))) (vector max_row max_val))) (def max_row2 (first ws1)) (def max_val2 (second ws1)) (if (>= row (dec (count a))) (vector max_row2 max_val2) (get_max_row (inc row) col a max_row2 max_val2))) ;ピボット選択 (defn pivoting [pivot a b a2 b2] (def ws1 (get_max_row 0 pivot a 0 0.0)) (def max_row (first ws1)) (def max_val (second ws1)) (def a3 (cons (nth a max_row) a2)) (def b3 (cons (nth b max_row) b2)) (def a4 (concat (take max_row a) (drop (inc max_row) a))) (def b4 (concat (take max_row b) (drop (inc max_row) b))) (if (>= pivot (dec N)) (vector (reverse a3) (reverse b3)) (pivoting (inc pivot) a4 b4 a3 b3))) ;前進消去 (defn forward_elim_loop [pivot row a b] (def s (/ (nth (nth a row) pivot) (nth (nth a pivot) pivot))) (def a_map (map (fn [a_piv a_row] (- a_row (* a_piv s))) (nth a pivot) (nth a row))) (def a2 (concat (take row a) (cons (vec a_map) (drop (inc row) a)))) (def x (- (nth b row) (* (nth b pivot) s))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (def next_row (if (= (inc row) pivot) (+ row 2) (inc row))) (if (< next_row N) (forward_elim_loop pivot next_row a2 b2) (vector a2 b2))) (defn forward_elimination [pivot a b] (def row (if (= pivot 0) 1 0)) (def ws1 (if (< pivot N) (forward_elim_loop pivot row a b) (vector a b))) (def a2 (first ws1)) (def b2 (second ws1)) (if (< pivot (dec N)) (forward_elimination (inc pivot) a2 b2) (vector a2 b2))) ;後退代入 (defn backward_substitution [row a b] (def x (/ (nth b row) (nth (nth a row) row))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (> row 0) (cons x (backward_substitution (dec row) a b2)) (cons x []))) ;ピボット選択 (def ws1 (pivoting 0 a b [] [])) (def a1 (first ws1)) (def b1 (vec (second ws1))) (println "pivoting") (println "A") (disp_matrix a1) (println "B") (disp_vector b1) (println "") ;前進消去 (def ws2 (forward_elimination 0 a1 b1)) (def a2 (first ws2)) (def b2 (vec (second ws2))) (println "forward elimination") (println "A") (disp_matrix a2) (println "B") (disp_vector b2) (println "") ;後退代入 (def x (backward_substitution (dec N) a2 b2)) (println "X") (disp_vector (reverse x))
Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj1004.clj pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000
Haskell
import Text.Printf import Debug.Trace import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 各列で 一番値が大きい行を 探す get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double) get_max_row row col a max_row max_val = let -- 一番値が大きい行 (max_row2, max_val2) = if max_val < abs(a!!row!!col) then (row, abs(a!!row!!col)) else (max_row, max_val) in if row >= length(a) - 1 then (max_row2, max_val2) else get_max_row (row+1) col a max_row2 max_val2 -- ピボット選択 pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double]) pivoting pivot a b a2 b2 = let (max_row, max_val) = get_max_row 0 pivot a 0 0.0 a3 = (a!!max_row):a2 b3 = (b!!max_row):b2 a4 = (take (max_row) a) ++ (drop (max_row+1) a) b4 = (take (max_row) b) ++ (drop (max_row+1) b) in if pivot >= (n - 1) then (reverse a3, reverse b3) else pivoting (pivot+1) a4 b4 a3 b3 -- 前進消去 forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elim_loop pivot row a b = let s = a!!row!!pivot / a!!pivot!!pivot a_zip = zip (a!!pivot) (a!!row) new_row = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip a2 = (take row a) ++ (new_row:(drop (row+1) a)) x = b!!row - b!!pivot * s b2 = (take row b) ++ (x:(drop (row+1) b)) next_row = if row + 1 == pivot then row + 2 else row + 1 in if next_row < n then forward_elim_loop pivot next_row a2 b2 else (a2, b2) forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elimination pivot a b = let row = if pivot == 0 then 1 else 0 (a2, b2) = if pivot < n then forward_elim_loop pivot row a b else (a, b) in if pivot < (n - 1) then forward_elimination (pivot+1) a2 b2 else (a2, b2) -- 後退代入 backward_substitution::Int->[[Double]]->[Double]->[Double] backward_substitution row a b = let x = b!!row / (a!!row!!row) b2 = (take row b) ++ (x:(drop (row+1) b)) in if row > 0 then x:(backward_substitution (row-1) a b2) else x:[] main = do let a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]] let b = [8,17,20,16::Double] -- ピボット選択 let (a1, b1) = pivoting 0 a b [] [] putStrLn "pivoting" putStrLn "A" disp_matrix a1 putStrLn "B" disp_vector b1 putStrLn "" -- 前進消去 let (a2, b2) = forward_elimination 0 a1 b1 putStrLn "forward elimination" putStrLn "A" disp_matrix a2 putStrLn "B" disp_vector b2 putStrLn "" -- 後退代入 let x = backward_substitution (n-1) a2 b2 putStrLn "X" disp_vector (reverse x)
Z:\>runghc Hs1004.hs pivoting A 9.0000000000 2.0000000000 1.0000000000 1.0000000000 2.0000000000 8.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 7.0000000000 -2.0000000000 1.0000000000 -1.0000000000 -2.0000000000 6.0000000000 B 20.0000000000 16.0000000000 8.0000000000 17.0000000000 forward elimination A 9.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 7.5555555556 0.0000000000 0.0000000000 0.0000000000 0.0000000000 6.5882352941 0.0000000000 0.0000000000 0.0000000000 0.0000000000 5.3750000000 B 9.0000000000 15.1111111111 19.7647058824 21.5000000000 X 1.0000000000 2.0000000000 3.0000000000 4.0000000000