さまざまな言語で数値計算
Only Do What Only You Can Do
逆反復法
逆反復法を使って, 最小固有値と, 対応する固有ベクトルを求める
正方行列 $ \boldsymbol{A} $ の逆行列を求め
- $ \boldsymbol{A}^{-1} \boldsymbol{x} = \displaystyle\frac{\boldsymbol{x}}{\lambda} $
逆行列を計算するより,
- $ \boldsymbol{Ax}_k = \boldsymbol{LUx}_k = \boldsymbol{x}_{k-1} $
逆べき乗法ともいう.
VBScript
Option Explicit Private Const N = 3 '逆ベキ乗法で最小固有値を求める Call Main Private Sub Main() Dim a: a = Array(_ Array(5.0, 4.0, 1.0, 1.0), _ Array(4.0, 5.0, 1.0, 1.0), _ Array(1.0, 1.0, 4.0, 2.0), _ Array(1.0, 1.0, 2.0, 4.0)) Dim x: x = Array(1.0 ,0.0 ,0.0 ,0.0) 'LU分解 Call forward_elimination(a) '逆ベキ乗法 Dim lambda: lambda = inverse(a, x) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" WScript.StdOut.WriteLine Right(Space(14) & FormatNumber(lambda, 10, -1, 0, 0), 14) WScript.StdOut.WriteLine "eigenvector" Call disp_vector(x) End Sub '逆ベキ乗法 Private Function inverse(ByVal a, ByRef x0) Dim i, j, k Dim lambda: lambda = 0.0 '正規化 (ベクトル x0 の長さを１にする) Call normarize(x0) Dim e0: e0 = 0.0 For i = 0 To N e0 = e0 + x0(i) Next For k = 1 To 200 '１次元配列を表示 WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_vector(x0) 'Ly = b から y を求める (前進代入) Dim y: y = Array(0.0, 0.0, 0.0, 0.0) Call forward_substitution(a,y,x0) 'Ux = y から x を求める (後退代入) Dim x1: x1 = Array(0.0, 0.0, 0.0, 0.0) Call backward_substitution(a,x1,y) '内積 Dim p0: p0 = 0.0 Dim p1: p1 = 0.0 For i = 0 To N p0 = p0 + x1(i) * x1(i) p1 = p1 + x1(i) * x0(i) Next '固有値 lambda = p1 / p0 '正規化 (ベクトル x1 の長さを１にする) Call normarize(x1) '収束判定 Dim e1: e1 = 0.0 For i = 0 To N e1 = e1 + x1(i) Next If Abs(e1 - e0) < 0.00000000001 Then Exit For x0 = x1 e0 = e1 Next inverse = lambda End Function 'LU分解 Private Sub forward_elimination(ByRef a) Dim pivot, row, col For pivot = 0 To (N - 1) For row = (pivot + 1) To N Dim s: s = a(row)(pivot) / a(pivot)(pivot) For col = pivot To N a(row)(col) = a(row)(col) - a(pivot)(col) * s 'これが 上三角行列 Next a(row)(pivot) = s 'これが 下三角行列 Next Next End Sub '前進代入 (Ly = b から y を求める) Private Sub forward_substitution(ByVal a, ByRef y, ByVal b) Dim row, col For row = 0 To N For col = 0 To row b(row) = b(row) - a(row)(col) * y(col) Next y(row) = b(row) Next End Sub '後退代入 (Ux = y から x を求める) Private Sub backward_substitution(ByVal a, ByRef x, ByVal y) Dim row, col For row = N To 0 Step -1 For col = N To (row + 1) Step - 1 y(row) = y(row) - a(row)(col) * x(col) Next x(row) = y(row) / a(row)(row) Next End Sub '１次元配列を表示 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 '正規化 (ベクトルの長さを１にする) Private Sub normarize(ByRef x()) Dim s: s = 0.0 Dim i For i = 0 To N s = s + x(i) * x(i) Next s = Sqr(s) For i = 0 To N x(i) = x(i) / s Next End Sub
Z:\>cscript //nologo Z:\1102.vbs 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
JScript
var N = 4 // 逆ベキ乗法で最小固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var x = [1.0 ,0.0 ,0.0 ,0.0] // LU分解 forward_elimination(a) // 逆ベキ乗法 var lambda = inverse(a, x) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") WScript.StdOut.WriteLine((" " + lambda.toFixed(10)).slice(-14)) WScript.StdOut.WriteLine("eigenvector") disp_vector(x) } // 逆ベキ乗法 function inverse(a, x0) { var lambda = 0.0 // 正規化 (ベクトル x0 の長さを１にする) normarize(x0) var e0 = 0.0 for (var i = 0; i < N; i++) e0 += x0[i] for (var k = 1; k <= 200; k++) { // １次元配列を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_vector(x0) // Ly=b から y を求める (前進代入) var b = [0,0,0,0] var y = [0,0,0,0] for (var i = 0; i < N; i++) b[i] = x0[i] forward_substitution(a,y,b) // Ux=y から x を求める (後退代入) var x1 = [0,0,0,0] backward_substitution(a,x1,y) // 内積 var p0 = 0.0 var p1 = 0.0 for (var i = 0; i < N; i++) { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p1 / p0 // 正規化 (ベクトル x1 の長さを１にする) normarize(x1) // 収束判定 var e1 = 0.0 for (var i = 0; i < N; i++) e1 += x1[i] if (Math.abs(e0 - e1) < 0.00000000001) break for (var i = 0; i < N; i++) x0[i] = x1[i] e0 = e1 } return lambda } // 前進消去 function forward_elimination(a) { for (var pivot = 0; pivot < N - 1; pivot++) { for (var row = pivot + 1; row < N; row++) { var s = a[row][pivot] / a[pivot][pivot] for (var col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s // これが 上三角行列 a[row][pivot] = s // これが 下三角行列 } } } // 前進代入 function forward_substitution(a, y, b) { for (var row = 0; row < N; row++) { for (var col = 0; col < row; col++) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // 後退代入 function backward_substitution(a, x, y) { for (var row = N - 1; row >= 0; row--) { for (var col = N - 1; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // １次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 正規化 (ベクトルの長さを１にする) function normarize(x) { var s = 0.0 for (var i = 0; i < N; i++) s += x[i] * x[i] s = Math.sqrt(s) for (var i = 0; i < N; i++) x[i] /= s }
Z:\>cscript //nologo Z:\1102.js 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
PowerShell
set-variable -name N -value 3 -option constant # LU分解 function forward_elimination($a) { foreach ($pivot in 0..($N - 1)) { foreach ($row in ($pivot + 1)..$N) { $s = $a[$row][$pivot] / $a[$pivot][$pivot] foreach ($col in $pivot..$N) { $a[$row][$col] -= $a[$pivot][$col] * $s # これが 上三角行列 } $a[$row][$pivot] = $s # これが 下三角行列 } } } # 前進代入 function forward_substitution($a, $y, $b) { foreach ($row in 0..$N) { foreach ($col in 0..$row) { $b[$row] -= $a[$row][$col] * $y[$col] } $y[$row] = $b[$row] } } # 後退代入 function backward_substitution($a, $x, $y) { foreach ($row in $N..0) { if ($row -lt $N) { foreach ($col in $N..($row + 1)) { $y[$row] -= $a[$row][$col] * $x[$col] } } $x[$row] = $y[$row] / $a[$row][$row] } } # １次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 正規化 (ベクトルの長さを１にする) function normarize($x) { $s = 0.0 foreach ($i in 0..$N) { $s += $x[$i] * $x[$i] } $s = [Math]::Sqrt($s) foreach ($i in 0..$N) { $x[$i] /= $s } } # 逆ベキ乗法 function inverse($a, $x0) { $lambda = 0.0 # 正規化 (ベクトル x0 の長さを１にする) normarize $x0 $e0 = 0.0 foreach ($i in 0..$N) { $e0 += $x0[$i] } foreach ($k in 1..200) { # １次元配列を表示 Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_vector $x0 # Ly = b から y を求める (前進代入) $y = (0.0, 0.0, 0.0, 0.0) $b = (0.0, 0.0, 0.0, 0.0) foreach ($i in 0..$N) { $b[$i] = $x0[$i] } forward_substitution $a $y $b # Ux = y から x を求める (後退代入) $x1 = (0.0, 0.0, 0.0, 0.0) backward_substitution $a $x1 $y # 内積 $p0 = 0.0 $p1 = 0.0 foreach ($i in 0..$N) { $p0 += $x1[$i] * $x1[$i] $p1 += $x1[$i] * $x0[$i] } # 固有値 $lambda = $p1 / $p0 # 正規化 (ベクトル x1 の長さを１にする) normarize $x1 # 収束判定 $e1 = 0.0 foreach ($i in 0..$N) { $e1 += $x1[$i] } if ([Math]::Abs($e0 - $e1) -lt 0.00000000001) { break } $lambda0 = $lambda1 foreach ($i in 0..$N) { $x0[$i] = $x1[$i] } $e0 = $e1 } $lambda } # 逆ベキ乗法で最小固有値を求める $a = (5.0, 4.0, 1.0, 1.0), (4.0, 5.0, 1.0, 1.0), (1.0, 1.0, 4.0, 2.0), (1.0, 1.0, 2.0, 4.0) $x = (1.0 ,0.0 ,0.0 ,0.0) # LU分解 forward_elimination $a # 逆ベキ乗法 $lambda = (inverse $a $x) Write-Host Write-Host "eigenvalue" Write-Host ([String]::Format("{0,14:F10}", $lambda)) Write-Host "eigenvector" disp_vector $x
Z:\>powershell -file Z:\1102.ps1 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 16 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Perl
use constant N => 3; # 逆ベキ乗法で最小固有値を求める main(); sub main { my @a = ([5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]); my @x = (1.0, 0.0, 0.0, 0.0); # LU分解 forward_elimination(\@a); # 逆ベキ乗法 Inverse iteration $lambda = inverse(\@a, \@x); print "\neigenvalue\n"; printf("%14.10f", $lambda); print "\neigenvector\n"; disp_vector(\@x); } # 逆ベキ乗法 Inverse iteration sub inverse { my ($a, $x0) = @_; my $lambda = 0.0; # 正規化 (ベクトル x0 の長さを１にする) normarize(\@$x0); my $e0 = 0.0; for $i (0..N) { $e0 += $$x0[$i]; } for $k (1..200) { # １次元配列を表示 printf("%3d\t", $k); disp_vector(\@$x0); # Ly = b から y を求める (前進代入) my @b = (0.0, 0.0, 0.0, 0.0); my @y; for $i (0..N) { $b[$i] = $$x0[$i]; } forward_substitution(\@$a, \@y, \@b); # Ux = y から x を求める (後退代入) my @x1; backward_substitution(\@$a, \@x1, \@y); # 内積 inner product my $p0 = 0.0; my $p1 = 0.0; for $i (0..N) { $p0 += $x1[$i] * $x1[$i]; $p1 += $x1[$i] * $$x0[$i]; } # 固有値 eigenvalue $lambda = $p1 / $p0; # 正規化 (ベクトル x1 の長さを１にする) normarize(\@x1); # 収束判定 my $e1 = 0.0; for $i (0..N) { $e1 += $x1[$i]; } last if (abs($e0 - $e1) < 0.00000000001); for $i (0..N) { $$x0[$i] = $x1[$i]; } $e0 = $e1; } $lambda; } # LU分解 sub forward_elimination { my ($a) = @_; for $pivot (0..(N - 1)) { for $row (($pivot + 1)..N) { my $s = $$a[$row][$pivot] / $$a[$pivot][$pivot]; for $col ($pivot..N) { $$a[$row][$col] -= $$a[$pivot][$col] * $s; # これが 上三角行列 } $$a[$row][$pivot] = $s; # これが 下三角行列 } } } # 前進代入 sub forward_substitution { my ($a, $y, $b) = @_; for $row (0..N) { for $col (0..$row) { $$b[$row] -= $$a[$row][$col] * $$y[$col]; } $$y[$row] = $$b[$row]; } } # 後退代入 sub backward_substitution { my ($a, $x, $y) = @_; for $row (reverse(0..N)) { for $col (reverse(($row + 1)..N)) { $$y[$row] -= $$a[$row][$col] * $$x[$col]; } $$x[$row] = $$y[$row] / $$a[$row][$row]; } } # １次元配列を表示 sub disp_vector { my ($x) = @_; for $i (0..N) { printf("%14.10f\t", $$x[$i]); } print "\n"; } # 正規化 (ベクトルの長さを１にする) sub normarize { my ($x) = @_; my $s = 0.0; for $i (0..N) { $s += $$x[$i] * $$x[$i]; } $s = sqrt($s); for $i (0..N) { $$x[$i] /= $s; } }
Z:\>perl Z:\1102.pl 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
PHP
<?php define("N", 3); # 逆ベキ乗法で最小固有値を求める main(); function main() { $a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]]; $x = [1.0, 0.0, 0.0, 0.0]; # LU分解 forward_elimination($a); # 逆ベキ乗法 $lambda = inverse($a, $x); print "\neigenvalue\n"; printf("%14.10f", $lambda); print "\neigenvector\n"; disp_vector($x); } # 逆ベキ乗法 function inverse(&$a, &$x0) { $lambda = 0.0; # 正規化 (ベクトル x0 の長さを１にする) normarize($x0); $e0 = 0.0; foreach (range(0, N) as $i) { $e0 += $x0[$i]; } foreach (range(1, 200) as $k) { # １次元配列を表示 printf("%3d\t", $k); disp_vector($x0); # Ly=b から y を求める (前進代入) $y = [0.0, 0.0, 0.0, 0.0]; forward_substitution($a, $y, $x0); # Ux=y から x を求める (後退代入) $x1 = [0.0, 0.0, 0.0, 0.0]; backward_substitution($a, $x1, $y); # 内積 $p0 = 0.0; $p1 = 0.0; foreach (range(0, N) as $i) { $p0 += $x1[$i] * $x1[$i]; $p1 += $x1[$i] * $x0[$i]; } # 固有値 $lambda = $p1 / $p0; # 正規化 (ベクトル x1 の長さを１にする) normarize($x1); # 収束判定 $e1 = 0.0; foreach (range(0, N) as $i) { $e1 += $x1[$i]; } if (abs($e0 - $e1) < 0.00000000001) break; foreach (range(0, N) as $i) { $x0[$i] = $x1[$i]; } $e0 = $e1; } return $lambda; } # LU分解 function forward_elimination(&$a) { foreach (range(0, N - 1) as $pivot) { foreach (range($pivot + 1, N) as $row) { $s = $a[$row][$pivot] / $a[$pivot][$pivot]; foreach (range($pivot, N) as $col) $a[$row][$col] -= $a[$pivot][$col] * $s; # これが 上三角行列 $a[$row][$pivot] = $s; # これが 下三角行列 } } } # 前進代入 function forward_substitution($a, &$y, $b) { foreach (range(0, N) as $row) { foreach (range(0, $row) as $col) $b[$row] -= $a[$row][$col] * $y[$col]; $y[$row] = $b[$row]; } } # 後退代入 function backward_substitution($a, &$x, $y) { foreach (range(N, 0) as $row) { if (($row + 1) <= N) foreach (range(N, $row + 1) as $col) $y[$row] -= $a[$row][$col] * $x[$col]; $x[$row] = $y[$row] / $a[$row][$row]; } } # １次元配列を表示 function disp_vector($row) { foreach ($row as $col) printf("%14.10f\t", $col); print "\n"; } # 正規化 (ベクトルの長さを１にする) function normarize(&$x) { $s = 0.0; foreach (range(0, N) as $i) { $s += $x[$i] * $x[$i]; } $s = sqrt($s); foreach (range(0, N) as $i) { $x[$i] /= $s; } } ?>
Z:\>php Z:\1102.php 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
Python
# coding: Shift_JIS import math N = 4 # 前進消去 def forward_elimination(a): for pivot in range(0, N - 1, 1): for row in range(pivot + 1, N, 1): s = a[row][pivot] / a[pivot][pivot] for col in range(pivot, N, 1): a[row][col] -= a[pivot][col] * s # これが 上三角行列 a[row][pivot] = s # これが 下三角行列 # 前進代入 def forward_substitution(a, y, b): for row in range(0, N, 1): for col in range(0, row, 1): b[row] -= a[row][col] * y[col] y[row] = b[row] # 後退代入 def backward_substitution(a, x, y): for row in range(N - 1, -1, -1): for col in range(N - 1, row, -1): y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] # １次元配列を表示 def disp_vector(row): for col in row: print "%14.10f\t" % col, print "" # 正規化 (ベクトルの長さを１にする) def normarize(x): s = 0.0 for i in range(0, N, 1): s += x[i] * x[i] s = math.sqrt(s) for i in range(0, N, 1): x[i] /= s # 逆ベキ乗法 def inverse(a, x0): lambda0 = 0.0 # 正規化 (ベクトル x0 の長さを１にする) normarize(x0) e0 = 0.0 for i in range(0, N, 1): e0 += x0[i] for k in range(1, 200, 1): # １次元配列を表示 print "%3d\t" % k, disp_vector(x0) # Ly=b から y を求める (前進代入) y = [0.0, 0.0, 0.0, 0.0] b = [0.0, 0.0, 0.0, 0.0] for i in range(0, N, 1): b[i] = x0[i] forward_substitution(a, y, b) # Ux=y から x を求める (後退代入) x1 = [ 0.0, 0.0, 0.0, 0.0] backward_substitution(a, x1, y) # 内積 p0 = 0.0 p1 = 0.0 for i in range(0, N, 1): p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] # 固有値 lambda0 = p1 / p0 # 正規化 (ベクトル x1 の長さを１にする) normarize(x1) # 収束判定 e1 = 0.0 for i in range(0, N, 1): e1 += x1[i] if (abs(e0 - e1) < 0.00000000001): break for i in range(0, N, 1): x0[i] = x1[i] e0 = e1 return lambda0 # 逆ベキ乗法で最小固有値を求める a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] x = [1.0 ,0.0 ,0.0 ,0.0] # LU分解 forward_elimination(a) # 逆ベキ乗法 lambda0 = inverse(a, x) print "" print "eigenvalue" print "%14.10f" % lambda0 print "eigenvector" disp_vector(x)
Z:\>python Z:\1102.py 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
Ruby
# coding: Shift_JIS N = 3 # 前進消去 def forward_elimination(a) 0.upto(N - 1) do |pivot| (pivot + 1).upto(N) do |row| s = a[row][pivot] / a[pivot][pivot] pivot.upto(N) do |col| a[row][col] -= a[pivot][col] * s # これが 上三角行列 end a[row][pivot] = s # これが 下三角行列 end end end # 前進代入 def forward_substitution(a, y, b) (0..N).each do |row| (0..row).each do |col| b[row] -= a[row][col] * y[col] end y[row] = b[row] end end # 後退代入 def backward_substitution(a, x, y) (0..N).reverse_each do |row| ((row + 1)..N).reverse_each do |col| y[row] -= a[row][col] * x[col] end x[row] = y[row] / a[row][row] end end # １次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 正規化 (ベクトルの長さを１にする) def normarize(x) s = 0.0 (0..N).each do |i| s += x[i] * x[i] end s = Math.sqrt(s) (0..N).each do |i| x[i] /= s end end # 逆ベキ乗法 def inverse(a, x0) lambda = 0.0 # 正規化 (ベクトル x0 の長さを１にする) normarize(x0) e0 = 0.0 (0..N).each do |i| e0 += x0[i] end (1..200).each do |k| # １次元配列を表示 printf("%3d\t", k) disp_vector(x0) # Ly=b から y を求める (前進代入) y = [0.0, 0.0, 0.0, 0.0] b = [0.0, 0.0, 0.0, 0.0] (0..N).each do |i| b[i] = x0[i] end forward_substitution(a, y, b) # Ux=y から x を求める (後退代入) x1 = [ 0.0, 0.0, 0.0, 0.0] backward_substitution(a, x1, y) # 内積 p0 = 0.0 p1 = 0.0 (0..N).each do |i| p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] end # 固有値 lambda = p1 / p0 # 正規化 (ベクトル x1 の長さを１にする) normarize(x1) # 収束判定 e1 = 0.0 (0..N).each do |i| e1 += x1[i] end break if ((e0 - e1).abs < 0.00000000001) (0..N).each do |i| x0[i] = x1[i] end e0 = e1 end lambda end # 逆ベキ乗法で最小固有値を求める a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] x = [1.0 ,0.0 ,0.0 ,0.0] # LU分解 forward_elimination(a) # 逆ベキ乗法 lambda0 = inverse(a, x) puts "" puts "eigenvalue" printf("%14.10f\n", lambda0) puts "eigenvector" disp_vector(x)
Z:\>ruby Z:\1102.rb 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
Groovy
N = 3 // 逆ベキ乗法で最小固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def x = [1.0 ,0.0 ,0.0 ,0.0] as double[] // LU分解 forward_elimination(a) // 逆ベキ乗法 lambda = inverse(a, x) println() println("eigenvalue") printf("%14.10f\n" , lambda) println("eigenvector") disp_vector(x) } // 逆ベキ乗法 def inverse(a, x0) { lambda = 0.0 // 正規化 (ベクトル x0 の長さを１にする) normarize(x0) e0 = 0.0 for (i in 0..N) e0 += x0[i] for (k in 1..200) { // １次元配列を表示 printf("%3d\t" , k) disp_vector(x0) // Ly=b から y を求める (前進代入) b = [0.0,0.0,0.0,0.0] y = [0.0,0.0,0.0,0.0] for (i in 0..N) b[i] = x0[i] forward_substitution(a,y,b) // Ux=y から x を求める (後退代入) x1 = [0.0,0.0,0.0,0.0] backward_substitution(a,x1,y) // 内積 p0 = 0.0 p1 = 0.0 for (i in 0..N) { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p1 / p0 // 正規化 (ベクトル x1 の長さを１にする) normarize(x1) // 収束判定 e1 = 0.0 for (i in 0..N) e1 += x1[i] if (Math.abs(e0 - e1) < 0.00000000001) break for (i in 0..N) x0[i] = x1[i] e0 = e1 } lambda } // 前進消去 def forward_elimination(a) { for (pivot in 0..(N - 1)) { for (row in (pivot + 1)..N) { s = a[row][pivot] / a[pivot][pivot] for (col in pivot..N) a[row][col] -= a[pivot][col] * s // これが 上三角行列 a[row][pivot] = s // これが 下三角行列 } } } // Ly=b から y を求める (前進代入) def forward_substitution(a, y, b) { for (row in 0..N) { for (col in 0..row) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // Ux=y から x を求める (後退代入) def backward_substitution(a, x, y) { for (row = N; row >= 0; row--) { for (col = N; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // １次元配列を表示 def disp_vector(row) { for (col in row) { printf("%14.10f\t" , col) } println() } // 正規化 (ベクトルの長さを１にする) def normarize(x) { s = 0.0 for (i in 0..N) s += x[i] * x[i] s = Math.sqrt(s) for (i in 0..N) x[i] /= s }
Z:\>gvy Gvy1102.gvy 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
Pascal
program Pas1102(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // １次元配列を表示 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; // 正規化 (ベクトルの長さを１にする) procedure normarize(var x:array of Double); var s:Double = 0.0; i:Integer; begin for i := Low(x) to High(x) do s := s + x[i] * x[i]; s := Sqrt(s); for i := Low(x) to High(x) do x[i] := x[i] / s; end; // LU分解 procedure forward_elimination(var a:TwoDimArray); var s:Double; pivot, row, col:Integer; begin for pivot := Low(a) to High(a) - 1 do begin for row := (pivot + 1) to High(a) do begin s := a[row,pivot] / a[pivot,pivot]; for col := pivot to High(a) do a[row,col] := a[row,col] - a[pivot,col] * s; // これが 上三角行列 a[row,pivot] := s; // これが 下三角行列 end; end; end; // 前進代入 (Ly = b から y を求める) procedure forward_substitution(a:TwoDimArray; var y:array of Double; b:array of Double); var row, col:Integer; begin for row := Low(a) to High(a) do begin for col := Low(a) to row do b[row] := b[row] - a[row,col] * y[col]; y[row] := b[row]; end; end; // 後退代入 (Ux = y から x を求める) procedure backward_substitution(a:TwoDimArray; var x:array of Double; var y:array of Double); var row, col:Integer; begin for row := High(a) downto Low(a) do begin for col := High(a) downto row + 1 do y[row] := y[row] - a[row,col] * x[col]; x[row] := y[row] / a[row,row]; end; end; // 逆ベキ乗法 function inverse(a:TwoDimArray; var x0:array of Double):Double; var lambda:Double = 0.0; p0, p1, e0, e1:Double; i, k:Integer; b:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); y:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); x1:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); begin // 正規化 (ベクトル x0 の長さを１にする) normarize(x0); e0 := 0.0; for i := Low(x0) to High(x0) do e0 := e0 + x0[i]; for k := 1 to 100 do begin // １次元配列を表示 write(format('%3d'#9, [k])); disp_vector(x0); // Ly = b から y を求める (前進代入) for i := Low(x0) to High(x0) do begin b[i] := x0[i]; y[i] := 0; end; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) for i := Low(x1) to High(x1) do x1[i] := 0; backward_substitution(a,x1,y); // 内積 p0 := 0.0; p1 := 0.0; for i := Low(x0) to High(x0) do begin p0 := p0 + x1[i] * x1[i]; p1 := p1 + x1[i] * x0[i]; end; // 固有値 lambda := p1 / p0; // 正規化 (ベクトル x1 の長さを１にする) normarize(x1); // 収束判定 e1 := 0.0; for i := Low(x0) to High(x0) do e1 := e1 + x1[i]; if Abs(e1 - e0) < 0.00000000001 then break; for i := Low(x0) to High(x0) do x0[i] := x1[i]; e0 := e1; end; inverse := lambda; end; // 逆ベキ乗法で最小固有値を求める var a :TwoDimArray = ((5.0, 4.0, 1.0, 1.0) ,(4.0, 5.0, 1.0, 1.0) ,(1.0, 1.0, 4.0, 2.0) ,(1.0, 1.0, 2.0, 4.0)); x :array [0..N] of Double = (1.0, 0.0, 0.0, 0.0); lambda:Double; begin // LU分解 forward_elimination(a); // 逆ベキ乗法 lambda := inverse(a, x); writeln(); writeln('eigenvalue'); writeln(format('%14.10f', [lambda])); writeln('eigenvector'); disp_vector(x); end.
Z:\>fpc -v0 -l- Pas1102.pp Z:\>Pas1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 16 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Ada
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1102 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 := ((5.0, 4.0, 1.0, 1.0) ,(4.0, 5.0, 1.0, 1.0) ,(1.0, 1.0, 4.0, 2.0) ,(1.0, 1.0, 2.0, 4.0)); x:Long_Float_Array := (1.0, 0.0, 0.0, 0.0); -- １次元配列を表示 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; -- 正規化 (ベクトルの長さを１にする) procedure normarize(x:in out Long_Float_Array) is s:Long_Float; begin s := 0.0; for i in x'Range loop s := s + x(i) * x(i); end loop; s := Sqrt(s); for i in x'Range loop x(i) := x(i) / s; end loop; end normarize; -- LU分解 procedure forward_elimination(a:in out Long_Float_TwoDimArray) is s:Long_Float; begin for pivot in a'First..(a'Last - 1) loop for row in (pivot + 1)..a'Last loop 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; a(row,pivot) := s; -- これが 下三角行列 end loop; end loop; end forward_elimination; -- 前進代入 (Ly = b から y を求める) procedure forward_substitution(a:Long_Float_TwoDimArray; y:in out Long_Float_Array; b:in out Long_Float_Array) is begin for row in a'Range loop for col in a'First..row loop b(row) := b(row) - a(row,col) * y(col); end loop; y(row) := b(row); end loop; end forward_substitution; -- 後退代入 (Ux = y から x を求める) procedure backward_substitution(a:Long_Float_TwoDimArray; x:in out Long_Float_Array; y:in out Long_Float_Array) is begin for row in reverse a'Range loop for col in reverse (row + 1)..a'Last loop y(row) := y(row) - a(row,col) * x(col); end loop; x(row) := y(row) / a(row,row); end loop; end backward_substitution; -- 逆ベキ乗法 procedure inverse(a:Long_Float_TwoDimArray; x0:in out Long_Float_Array; lambda:out Long_Float) is p0, p1, e0, e1:Long_Float; b:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); y:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); x1:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); begin -- 正規化 (ベクトル x0 の長さを１にする) normarize(x0); e0 := 0.0; for i in x0'Range loop e0 := e0 + x0(i); end loop; for k in 1..100 loop -- １次元配列を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_vector(x0); -- Ly = b から y を求める (前進代入) for i in x0'Range loop b(i) := x0(i); y(i) := 0.0; end loop; forward_substitution(a,y,b); -- Ux = y から x を求める (後退代入) for i in x1'Range loop x1(i) := 0.0; end loop; backward_substitution(a,x1,y); -- 内積 p0 := 0.0; p1 := 0.0; for i in x0'Range loop p0 := p0 + x1(i) * x1(i); p1 := p1 + x1(i) * x0(i); end loop; -- 固有値 lambda := p1 / p0; -- 正規化 (ベクトル x1 の長さを１にする) normarize(x1); -- 収束判定 e1 := 0.0; for i in x0'Range loop e1 := e1 + x1(i); end loop; if Abs(e1 - e0) < 0.00000000001 then exit; end if; for i in x0'Range loop x0(i) := x1(i); end loop; e0 := e1; end loop; end inverse; -- 逆ベキ乗法で最小固有値を求める lambda:Long_Float := 0.0; begin -- LU分解 forward_elimination(a); -- 逆ベキ乗法 inverse(a, x, lambda); Put_Line(""); Put_Line("eigenvalue"); Put(lambda, Fore=>3, Aft=>10, Exp=>0); New_Line; Put_Line("eigenvector"); disp_vector(x); end Ada1102;
xxxxxx@yyyyyy /Z $ gnatmake Ada1102.adb xxxxxx@yyyyyy /Z $ Ada1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
VB.NET
Option Explicit Module VB1102 Private Const N As Integer = 3 '逆ベキ乗法で最小固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} Dim x() As Double = {1.0 ,0.0 ,0.0 ,0.0} 'LU分解 forward_elimination(a) '逆ベキ乗法 Dim lambda As Double = inverse(a, x) Console.WriteLine() Console.WriteLine("eigenvalue") Console.WriteLine(String.Format("{0,14:F10}", lambda)) Console.WriteLine("eigenvector") disp_vector(x) End Sub '逆ベキ乗法 Private Function inverse(ByVal a(,) As Double, ByRef x0() As Double) As Double Dim lambda As Double = 0.0 '正規化 (ベクトル x0 の長さを１にする) normarize(x0) Dim e0 As Double = 0.0 For i As Integer = 0 To N e0 += x0(i) Next For k As Integer = 1 To 100 '１次元配列を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_vector(x0) 'Ly = b から y を求める (前進代入) Dim b() As Double = New Double(N){} Dim y() As Double = New Double(N){} For i As Integer = 0 To N b(i) = x0(i) Next forward_substitution(a,y,b) 'Ux = y から x を求める (後退代入) Dim x1() As Double = New Double(N){} backward_substitution(a,x1,y) '内積 Dim p0 As Double = 0.0 Dim p1 As Double = 0.0 For i As Integer = 0 To N p0 += x1(i) * x1(i) p1 += x1(i) * x0(i) Next '固有値 lambda = p1 / p0 '正規化 (ベクトル x1 の長さを１にする) normarize(x1) '収束判定 Dim e1 As Double = 0.0 For i As Integer = 0 To N e1 += x1(i) Next If Math.Abs(e1 - e0) < 0.00000000001 Then Exit For For i As Integer = 0 To N x0(i) = x1(i) Next e0 = e1 Next Return lambda End Function 'LU分解 Private Sub forward_elimination(ByRef a(,) As Double) For pivot As Integer = 0 To (N - 1) For row As Integer = (pivot + 1) To N 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 a(row, pivot) = s 'これが 下三角行列 Next Next End Sub '前進代入 (Ly = b から y を求める) Private Sub forward_substitution(ByVal a(,) As Double, ByRef y() As Double, ByVal b() As Double) For row As Integer = 0 To N For col As Integer = 0 To row b(row) -= a(row, col) * y(col) Next y(row) = b(row) Next End Sub '後退代入 (Ux = y から x を求める) Private Sub backward_substitution(ByVal a(,) As Double, ByRef x() As Double, ByVal y() As Double) For row As Integer = N To 0 Step -1 For col As Integer = N To (row + 1) Step - 1 y(row) -= a(row, col) * x(col) Next x(row) = y(row) / a(row, row) Next End Sub '１次元配列を表示 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 '正規化 (ベクトルの長さを１にする) Private Sub normarize(ByRef x() As Double) Dim s As Double = 0.0 For i As Integer = 0 To N s += x(i) * x(i) Next s = Math.Sqrt(s) For i As Integer = 0 To N x(i) /= s Next End Sub End Module
Z:\>vbc -nologo VB1102.vb Z:\>VB1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 16 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
C#
using System; public class CS1102 { private const int N = 4; // 逆ベキ乗法で最小固有値を求める public static void Main() { double[,] a = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double[] x = {1.0 ,0.0 ,0.0 ,0.0}; // LU分解 forward_elimination(a); // 逆ベキ乗法 double lambda = inverse(a, x); Console.WriteLine(); Console.WriteLine("eigenvalue"); Console.WriteLine(string.Format("{0,14:F10}", lambda)); Console.WriteLine("eigenvector"); disp_vector(x); } // 逆ベキ乗法 private static double inverse(double[,] a, double[] x0) { double lambda = 0.0; // 正規化 (ベクトル x0 の長さを１にする) normarize(x0); double e0 = 0.0; for (int i = 0; i < N; i++) e0 += x0[i]; for (int k = 1; k <= 200; k++) { // １次元配列を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_vector(x0); // Ly = b から y を求める (前進代入) double[] b = new double[N]; double[] y = new double[N]; for (int i = 0; i < N; i++) b[i] = x0[i]; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) double[] x1 = new double[N]; backward_substitution(a,x1,y); // 内積 double p0 = 0.0; double p1 = 0.0; for (int i = 0; i < N; i++) { p0 += x1[i] * x1[i]; p1 += x1[i] * x0[i]; } // 固有値 lambda = p1 / p0; // 正規化 (ベクトル x1 の長さを１にする) normarize(x1); // 収束判定 double e1 = 0.0; for (int i = 0; i < N; i++) e1 += x1[i]; if (Math.Abs(e0 - e1) < 0.00000000001) break; for (int i = 0; i < N; i++) x0[i] = x1[i]; e0 = e1; } return lambda; } // LU分解 private static void forward_elimination(double[,] a) { for (int pivot = 0; pivot < N - 1; pivot++) { for (int row = pivot + 1; row < N; row++) { double s = a[row,pivot] / a[pivot,pivot]; for (int col = pivot; col < N; col++) a[row,col] -= a[pivot,col] * s; // これが 上三角行列 a[row,pivot] = s; // これが 下三角行列 } } } // 前進代入 private static void forward_substitution(double[,] a, double[] y, double[] b) { for (int row = 0; row < N; row++) { for (int col = 0; col < row; col++) b[row] -= a[row,col] * y[col]; y[row] = b[row]; } } // 後退代入 private static void backward_substitution(double[,] a, double[] x, double[] y) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) y[row] -= a[row,col] * x[col]; x[row] = y[row] / a[row,row]; } } // １次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 正規化 (ベクトルの長さを１にする) private static void normarize(double[] x) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = Math.Sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; } }
Z:\>csc -nologo CS1102.cs Z:\>CS1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 16 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Java
import java.lang.*; public class Java1102 { private static final int N = 4; // 逆ベキ乗法で最小固有値を求める public static void main(String []args) { double[][] a = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double[] x = {1.0 ,0.0 ,0.0 ,0.0}; // LU分解 forward_elimination(a); // 逆ベキ乗法 double lambda = inverse(a, x); System.out.println(); System.out.println("eigenvalue"); System.out.println(String.format("%14.10f", lambda)); System.out.println("eigenvector"); disp_vector(x); } // 逆ベキ乗法 private static double inverse(double[][] a, double[] x0) { double lambda = 0.0; // 正規化 (ベクトル x0 の長さを１にする) normarize(x0); double e0 = 0.0; for (int i = 0; i < N; i++) e0 += x0[i]; for (int k = 1; k <= 100; k++) { // 正規化 (ベクトル x0 の長さを１にする) normarize(x0); // １次元配列を表示 System.out.print(String.format("%3d\t", k)); disp_vector(x0); // Ly = b から y を求める (前進代入) double[] b = new double[N]; double[] y = new double[N]; for (int i = 0; i < N; i++) b[i] = x0[i]; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) double[] x1 = new double[N]; backward_substitution(a,x1,y); // 内積 double p0 = 0.0; double p1 = 0.0; for (int i = 0; i < N; i++) { p0 += x1[i] * x1[i]; p1 += x1[i] * x0[i]; } // 固有値 lambda = p1 / p0; // 正規化 (ベクトル x1 の長さを１にする) normarize(x1); // 収束判定 double e1 = 0.0; for (int i = 0; i < N; i++) e1 += x1[i]; if (Math.abs(e0 - e1) < 0.00000000001) break; for (int i = 0; i < N; i++) x0[i] = x1[i]; e0 = e1; } return lambda; } // LU分解 private static void forward_elimination(double[][] a) { for (int pivot = 0; pivot < N - 1; pivot++) { for (int row = pivot + 1; row < N; row++) { double s = a[row][pivot] / a[pivot][pivot]; for (int col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s; // これが 上三角行列 a[row][pivot] = s; // これが 下三角行列 } } } // 前進代入 private static void forward_substitution(double[][] a, double[] y, double[] b) { for (int row = 0; row < N; row++) { for (int col = 0; col < row; col++) b[row] -= a[row][col] * y[col]; y[row] = b[row]; } } // 後退代入 private static void backward_substitution(double[][] a, double[] x, double[] y) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) y[row] -= a[row][col] * x[col]; x[row] = y[row] / a[row][row]; } } // １次元配列を表示 private static void disp_vector(double[] row) { for (double col: row) System.out.print(String.format("%14.10f\t", col)); System.out.println(); } // 正規化 (ベクトルの長さを１にする) private static void normarize(double[] x) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = Math.sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; } }
Z:\>javac Java1102.java Z:\>java Java1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
C++
#include <iostream> #include <iomanip> #include <math.h> using namespace std; const int N = 4; // 逆ベキ乗法 double inverse(double a[N][N], double x0[N]); // LU分解 void forward_elimination(double a[N][N]); // 前進代入 void forward_substitution(double a[N][N], double y[N], double b[N]); // 後退代入 void backward_substitution(double a[N][N], double x[N], double y[N]); // １次元配列を表示 void disp_vector(double row[N]); // 正規化 (ベクトルの長さを１にする) void normarize(double x[N]); // 逆ベキ乗法で最小固有値を求める int main() { double a[N][N] = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double x[N] = {1.0 ,0.0 ,0.0 ,0.0}; // LU分解 forward_elimination(a); // 逆ベキ乗法 double lambda = inverse(a, x); cout << endl << "eigenvalue" << endl; cout << setw(14) << fixed << setprecision(10) << lambda; cout << endl << "eigenvector" << endl; disp_vector(x); return 0; } // 逆ベキ乗法 double inverse(double a[N][N], double x0[N]) { double lambda = 0.0; // 正規化 (ベクトル x0 の長さを１にする) normarize(x0); double e0 = 0.0; for (int i = 0; i < N; i++) e0 += x0[i]; for (int k = 1; k <= 100; k++) { // １次元配列を表示 cout << setw(3) << k << "\t"; disp_vector(x0); // Ly = b から y を求める (前進代入) double b[N] = {0.0 ,0.0 ,0.0 ,0.0}; double y[N] = {0.0 ,0.0 ,0.0 ,0.0}; for (int i = 0; i < N; i++) b[i] = x0[i]; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) double x1[N] = {0.0 ,0.0 ,0.0 ,0.0}; backward_substitution(a,x1,y); // 内積 double p0 = 0.0; double p1 = 0.0; for (int i = 0; i < N; i++) { p0 += x1[i] * x1[i]; p1 += x1[i] * x0[i]; } // 固有値 lambda = p1 / p0; // 正規化 (ベクトル x1 の長さを１にする) normarize(x1); // 収束判定 double e1 = 0.0; for (int i = 0; i < N; i++) e1 += x1[i]; if (fabs(e0 - e1) < 0.00000000001) break; for (int i = 0; i < N; i++) x0[i] = x1[i]; e0 = e1; } return lambda; } // LU分解 void forward_elimination(double a[N][N]) { for (int pivot = 0; pivot < N - 1; pivot++) { for (int row = pivot + 1; row < N; row++) { double s = a[row][pivot] / a[pivot][pivot]; for (int col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s; // これが 上三角行列 a[row][pivot] = s; // これが 下三角行列 } } } // 前進代入 void forward_substitution(double a[N][N], double y[N], double b[N]) { for (int row = 0; row < N; row++) { for (int col = 0; col < row; col++) b[row] -= a[row][col] * y[col]; y[row] = b[row]; } } // 後退代入 void backward_substitution(double a[N][N], double x[N], double y[N]) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) y[row] -= a[row][col] * x[col]; x[row] = y[row] / a[row][row]; } } // １次元配列を表示 void disp_vector(double row[N]) { for (int i = 0; i < N; i++) cout << setw(14) << fixed << setprecision(10) << row[i] << "\t"; cout << endl; } // 正規化 (ベクトルの長さを１にする) void normarize(double x[N]) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; }
Z:\>bcc32 -q CP1102.cpp cp1102.cpp: Z:\>CP1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
Objective-C
#import <Foundation/Foundation.h> #import <math.h> #define N 4 // 逆ベキ乗法 double inverse(double a[N][N], double x0[N]); // LU分解 void forward_elimination(double a[N][N]); // 前進代入 void forward_substitution(double a[N][N], double y[N], double b[N]); // 後退代入 void backward_substitution(double a[N][N], double x[N], double y[N]); // １次元配列を表示 void disp_vector(double row[N]); // 正規化 (ベクトルの長さを１にする) void normarize(double x[N]); // 逆ベキ乗法で最小固有値を求める int main() { double a[N][N] = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double x[N] = {1.0 ,0.0 ,0.0 ,0.0}; // LU分解 forward_elimination(a); // 逆ベキ乗法 Inverse iteration double lambda = inverse(a, x); printf("\neigenvalue\n"); printf("%14.10f\n", lambda); printf("eigenvector\n"); disp_vector(x); return 0; } // 逆ベキ乗法 double inverse(double a[N][N], double x0[N]) { int i, j, k; double lambda = 0.0; // 正規化 (ベクトル x0 の長さを１にする) normarize(x0); double e0 = 0.0; for (i = 0; i < N; i++) e0 += x0[i]; for (k = 1; k <= 200; k++) { // １次元配列を表示 printf("%3d\t", k); disp_vector(x0); // Ly = b から y を求める (前進代入) double b[N] = {0.0 ,0.0 ,0.0 ,0.0}; double y[N] = {0.0 ,0.0 ,0.0 ,0.0}; for (i = 0; i < N; i++) b[i] = x0[i]; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) double x1[N] = {0.0 ,0.0 ,0.0 ,0.0}; backward_substitution(a,x1,y); // 内積 double p0 = 0.0; double p1 = 0.0; for (i = 0; i < N; i++) { p0 += x1[i] * x1[i]; p1 += x1[i] * x0[i]; } // 固有値 lambda = p1 / p0; // 正規化 (ベクトル x1 の長さを１にする) normarize(x1); double e1 = 0.0; // 収束判定 for (i = 0; i < N; i++) e1 += x1[i]; if (fabs(e0 - e1) < 0.00000000001) break; for (i = 0; i < N; i++) x0[i] = x1[i]; e0 = e1; } return lambda; } // LU分解 void forward_elimination(double a[N][N]) { int pivot, row, col; for (pivot = 0; pivot < N - 1; pivot++) { for (row = pivot + 1; row < N; row++) { double s = a[row][pivot] / a[pivot][pivot]; for (col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s; // これが 上三角行列 a[row][pivot] = s; // これが 下三角行列 } } } // 前進代入 void forward_substitution(double a[N][N], double y[N], double b[N]) { int row, col; for (row = 0; row < N; row++) { for (col = 0; col < row; col++) b[row] -= a[row][col] * y[col]; y[row] = b[row]; } } // 後退代入 void backward_substitution(double a[N][N], double x[N], double y[N]) { int row, col; for (row = N - 1; row >= 0; row--) { for (col = N - 1; col > row; col--) y[row] -= a[row][col] * x[col]; x[row] = y[row] / a[row][row]; } } // １次元配列を表示 void disp_vector(double row[N]) { int i; for (i = 0; i < N; i++) printf("%14.10f\t", row[i]); printf("\n"); } // 正規化 (ベクトルの長さを１にする) void normarize(double x[N]) { int i; double s = 0.0; for (i = 0; i < N; i++) s += x[i] * x[i]; s = sqrt(s); for (i = 0; i < N; i++) x[i] /= s; }
xxxxxx@yyyyyy /Z $ gcc -o OC1102 OC1102.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
D
import std.stdio; import std.math; const int N = 4; // 逆ベキ乗法で最小固有値を求める void main(string[] args) { double a[N][N] = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]]; double x[N] = [1.0 ,0.0 ,0.0 ,0.0]; // LU分解 forward_elimination(a); // 逆ベキ乗法 double lambda = inverse(a, x); writefln("\neigenvalue"); writefln("%14.10f", lambda); writefln("eigenvector"); disp_vector(x); } // 逆ベキ乗法 double inverse(double[N][N] a, ref double[N] x0) { double lambda = 0.0; // 正規化 (ベクトル x0 の長さを１にする) normarize(x0); double e0 = 0.0; foreach (i; 0..N) e0 += x0[i]; foreach (k; 1..200) { // １次元配列を表示 writef("%3d\t", k); disp_vector(x0); // Ly = b から y を求める (前進代入) double[N] b = [0.0 ,0.0 ,0.0 ,0.0]; double[N] y = [0.0 ,0.0 ,0.0 ,0.0]; foreach (i; 0..N) b[i] = x0[i]; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) double[N] x1 = [0.0 ,0.0 ,0.0 ,0.0]; backward_substitution(a,x1,y); // 内積 double p0 = 0.0; double p1 = 0.0; foreach (i; 0..N) { p0 += x1[i] * x1[i]; p1 += x1[i] * x0[i]; } // 固有値 lambda = p1 / p0; // 正規化 (ベクトル x1 の長さを１にする) normarize(x1); // 収束判定 double e1 = 0.0; foreach (i; 0..N) e1 += x1[i]; if (fabs(e0 - e1) < 0.00000000001) break; foreach (i; 0..N) x0[i] = x1[i]; e0 = e1; } return lambda; } // LU分解 void forward_elimination(ref double[N][N] a) { foreach (pivot; 0..N) { foreach (row; (pivot + 1)..N) { double s = a[row][pivot] / a[pivot][pivot]; foreach (col; pivot..N) a[row][col] -= a[pivot][col] * s; // これが 上三角行列 a[row][pivot] = s; // これが 下三角行列 } } } // 前進代入 void forward_substitution(double[N][N] a, ref double[N] y, ref double[N] b) { foreach (row; 0..N) { foreach (col; 0..row) b[row] -= a[row][col] * y[col]; y[row] = b[row]; } } // 後退代入 void backward_substitution(double[N][N] a, ref double[N] x, ref double[N] y) { foreach_reverse (row; 0..N) { foreach_reverse (col; (row+1)..N) y[row] -= a[row][col] * x[col]; x[row] = y[row] / a[row][row]; } } // １次元配列を表示 void disp_vector(double[N] row) { foreach(col; row) writef("%14.10f\t", col); writefln(""); } // 正規化 (ベクトルの長さを１にする) void normarize(ref double[N] x) { double s = 0.0; foreach (i; 0..N) s += x[i] * x[i]; s = sqrt(s); foreach (i; 0..N) x[i] /= s; }
Z:\>dmd D1102.d Z:\>D1102 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606766 -0.0280606766 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242639 -0.0004242639 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178190 -0.0000178190 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000057 -0.0000000057 13 0.7071067818 -0.7071067806 -0.0000000011 -0.0000000011 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
Go
package main import "fmt" import "math" const N = 4 func main() { var a [N][N]float64 = [N][N]float64{{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} var x = []float64 {1.0 ,0.0 ,0.0 ,0.0} // LU分解 forward_elimination(&a) // 逆ベキ乗法 var lambda float64 = inverse(a, x) fmt.Println("\neigenvalue") fmt.Printf("%14.10f\n", lambda) fmt.Println("eigenvector") disp_vector(x) } // 逆ベキ乗法 func inverse(a[N][N]float64, x0[]float64) float64 { var lambda float64 = 0.0 // 正規化 (ベクトル x0 の長さを１にする) normarize(x0) var e0 float64 = 0.0 for i := 0; i < N; i++ { e0 += x0[i] } for k := 1; k < 200; k++ { // １次元配列を表示 fmt.Printf("%3d\t", k) disp_vector(x0) // Ly = b から y を求める (前進代入) var b = []float64 {0.0 ,0.0 ,0.0 ,0.0} var y = []float64 {0.0 ,0.0 ,0.0 ,0.0} for i := 0; i < N; i++ { b[i] = x0[i] } forward_substitution(a,y,b) // Ux = y から x を求める (後退代入) var x1 = []float64 {0.0 ,0.0 ,0.0 ,0.0} backward_substitution(a,x1,y) // 内積 var p0 float64 = 0.0 var p1 float64 = 0.0 for i := 0; i < N; i++ { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p1 / p0 // 正規化 (ベクトル x1 の長さを１にする) normarize(x1) // 収束判定 var e1 float64 = 0.0 for i := 0; i < N; i++ { e1 += x1[i] } if (math.Abs(e0 - e1) < 0.00000000001) { break } for i := 0; i < N; i++ { x0[i] = x1[i] } e0 = e1 } return lambda } // LU分解 func forward_elimination(a *[N][N]float64) { for pivot := 0; pivot < N - 1; pivot++ { for row := pivot + 1; row < N; row++ { var s = a[row][pivot] / a[pivot][pivot] for col := pivot; col < N; col++ { a[row][col] -= a[pivot][col] * s // これが 上三角行列 } a[row][pivot] = s // これが 下三角行列 } } } // 前進代入 func forward_substitution(a [N][N]float64, y []float64, b []float64) { for row := 0; row < N; row++ { for col := 0; col < row; col++ { b[row] -= a[row][col] * y[col] } y[row] = b[row] } } // 後退代入 func backward_substitution(a [N][N]float64, x []float64, y []float64) { for row := N - 1; row >= 0; row-- { for col := N - 1; col > row; col-- { y[row] -= a[row][col] * x[col] } x[row] = y[row] / a[row][row] } } // １次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 正規化 (ベクトルの長さを１にする) func normarize(x[]float64) { var s float64 = 0.0 for i := 0; i < N; i++ { s += x[i] * x[i] } s = math.Sqrt(s) for i := 0; i < N; i++ { x[i] /= s } }
Z:\>go run GO1102.go 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
Scala
object Scala1102 { val N = 3 // 逆ベキ乗法で最小固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) var x:Array[Double] = Array(1.0 ,0.0 ,0.0 ,0.0) // LU分解 forward_elimination(a) // 逆ベキ乗法 val (lambda, x1) = inverse(a, x) println() println("eigenvalue") println("%14.10f".format(lambda)) println("eigenvector") disp_vector(x1) } // 逆ベキ乗法 def inverse(a:Array[Array[Double]], x:Array[Double]): (Double, Array[Double]) = { var x0:Array[Double] = normarize(x) var x1:Array[Double] = Array(0.0 ,0.0 ,0.0 ,0.0) var lambda:Double = 0.0 var e:Double = 0.0 var k = 1 do { e = x0.sum // １次元配列を表示 print("%3d\t".format(k)) disp_vector(x0) k += 1 // Ly = b から y を求める (前進代入) var b = x0.clone var y:Array[Double] = Array(0.0 ,0.0 ,0.0 ,0.0) forward_substitution(a,y,b) // Ux = y から x を求める (後退代入) x1 = Array(0.0 ,0.0 ,0.0 ,0.0) backward_substitution(a,x1,y) // 内積 var (t1, t2) = (x0 zip x1).map(t => t match {case (f, s) => (s * s, s * f)}).unzip // 固有値 lambda = t2.sum / t1.sum // 正規化 (ベクトル x1 の長さを１にする) x1 = normarize(x1) x0 = x1.clone } while (k <= 100 && Math.abs(e - x1.sum) >= 0.00000000001) return (lambda, x0) } // LU分解 def forward_elimination(a:Array[Array[Double]]) { for (pivot <- 0 to N - 1) { for (row <- pivot + 1 to N) { var s:Double = a(row)(pivot) / a(pivot)(pivot) for (col <- pivot to N) a(row)(col) -= a(pivot)(col) * s // これが 上三角行列 a(row)(pivot) = s // これが 下三角行列 } } } // 前進代入 def forward_substitution(a:Array[Array[Double]], y:Array[Double], b:Array[Double]) { for (row <- 0 to N) { for (col <- 0 to row - 1) b(row) -= a(row)(col) * y(col) y(row) = b(row) } } // 後退代入 def backward_substitution(a:Array[Array[Double]], x:Array[Double], y:Array[Double]) { for (row <- (0 to N).reverse) { for (col <- (row + 1 to N).reverse) y(row) -= a(row)(col) * x(col) x(row) = y(row) / a(row)(row) } } // １次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 正規化 (ベクトルの長さを１にする) def normarize(x:Array[Double]):Array[Double] = { var s:Double = x.map(n => n * n).sum s = Math.sqrt(s) x.map(n => n / s) } }
Z:\>scala Scala1102.scala 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 16 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 -0.0000000000 -0.0000000000
F#
module Fs1102 open System let N = 3 // LU分解 let forward_elimination (a:float[][]) = for pivot in [0..N - 1] do for row in [pivot + 1 .. N] do 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 // これが 上三角行列 a.[row].[pivot] <- s // これが 下三角行列 // 前進代入 let forward_substitution (a:float[][]) (y:float[]) (b:float[]) = for row in [0..N] do for col in 0..(row - 1) do b.[row] <- b.[row] - a.[row].[col] * y.[col] y.[row] <- b.[row] // 後退代入 let backward_substitution (a:float[][]) (x:float[]) (y:float[]) = for row in (List.rev [0..N]) do for col in (List.rev [row + 1 .. N]) do y.[row] <- y.[row] - a.[row].[col] * x.[col] x.[row] <- y.[row] / a.[row].[row] // １次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 正規化 (ベクトルの長さを１にする) let normarize(x:float[]):float[] = let s = x |> Array.map (fun n -> n * n) |> Array.sum x |> Array.map(fun n -> n / Math.Sqrt(s)) // 逆ベキ乗法 let inverse (a:float[][]) (x:float[]) : float * float[] = let mutable x0:float[] = normarize x let mutable x1:float[] = [|0.0 ;0.0 ;0.0 ;0.0|] let mutable lambda = 0.0 let mutable e = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 100 && Math.Abs(e - (x1 |> Array.sum)) >= 0.00000000001) do if sw then sw <- false e <- x0 |> Array.sum // １次元配列を表示 printf "%3d\t" k disp_vector x0 k <- k + 1 // Ly = b から y を求める (前進代入) let mutable b = Array.copy x0 let mutable y:float[] = [|0.0 ;0.0 ;0.0 ;0.0|] forward_substitution a y b // Ux = y から x を求める (後退代入) let mutable x1 = [|0.0 ;0.0 ;0.0 ;0.0|] backward_substitution a x1 y // 内積 let t1, t2 = Array.zip x0 x1 |> Array.map(fun t -> match t with | (f, s) -> (s * s, s * f)) |> Array.unzip // 固有値 lambda <- (t2 |> Array.sum) / (t1 |> Array.sum) // 正規化 (ベクトル x1 の長さを１にする) x1 <- normarize x1 x0 <- Array.copy x1 (lambda, x0) // 逆ベキ乗法で最小固有値を求める let mutable a:float[][] = [| [|5.0; 4.0; 1.0; 1.0|]; [|4.0; 5.0; 1.0; 1.0|]; [|1.0; 1.0; 4.0; 2.0|]; [|1.0; 1.0; 2.0; 4.0|] |] let mutable x:float[] = [|1.0 ;0.0 ;0.0 ;0.0|] // LU分解 forward_elimination a // 逆ベキ乗法 let (lambda, x1) = inverse a x printfn "" printfn "eigenvalue" printfn "%14.10f" lambda printfn "eigenvector" disp_vector x1 exit 0
Z:\>fsi --nologo --quiet Fs1102.fs 1 1.0000000000 0.0000000000 0.0000000000 0.0000000000 2 0.7856989466 -0.6173348866 -0.0280606767 -0.0280606767 3 0.7182768487 -0.6956539558 -0.0084835848 -0.0084835848 4 0.7087990176 -0.7054049282 -0.0019798855 -0.0019798855 5 0.7073894400 -0.7068237547 -0.0004242640 -0.0004242640 6 0.7071576856 -0.7070558622 -0.0000876812 -0.0000876812 7 0.7071163975 -0.7070971642 -0.0000178191 -0.0000178191 8 0.7071086479 -0.7071049144 -0.0000035921 -0.0000035921 9 0.7071071489 -0.7071064135 -0.0000007212 -0.0000007212 10 0.7071068542 -0.7071067082 -0.0000001445 -0.0000001445 11 0.7071067957 -0.7071067666 -0.0000000289 -0.0000000289 12 0.7071067841 -0.7071067783 -0.0000000058 -0.0000000058 13 0.7071067818 -0.7071067806 -0.0000000012 -0.0000000012 14 0.7071067813 -0.7071067811 -0.0000000002 -0.0000000002 15 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 16 0.7071067812 -0.7071067812 0.0000000000 0.0000000000 eigenvalue 1.0000000000 eigenvector 0.7071067812 -0.7071067812 0.0000000000 0.0000000000