さまざまな言語で数値計算
Only Do What Only You Can Do
LR分解法
LR分解を使って, 全ての固有値を求める
対称行列 $ \boldsymbol{A} $ を 左三角行列 $ \boldsymbol{L} $ と 右三角行列 $ \boldsymbol{R} $ の積に分解し,
- $ \boldsymbol{A}_k = \boldsymbol{L}_k \boldsymbol{R}_k $
- $ \boldsymbol{A}_{k+1} = \boldsymbol{R}_k \boldsymbol{L}_k $
VBScript
Option Explicit Private Const N = 3 'LR分解で固有値を求める 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 l: l = Array(_ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0)) Dim u: u = Array(_ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0)) Dim i, j, k For k = 1 To 200 'LU分解 Call decomp(a, l, u) '行列の積 Call multiply(u, l, a) '対角要素を表示 WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_eigenvalue(a) '収束判定 Dim e: e = 0.0 For i = 1 To N For j = 0 To i - 1 e = e + Abs(a(i)(j)) Next Next If e < 0.00000000001 Then Exit For Next WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" Call disp_eigenvalue(a) End Sub '対角要素を表示 Private Sub disp_eigenvalue(ByVal a) Dim i For i = 0 To N WScript.StdOut.Write Right(Space(14) & FormatNumber(a(i)(i), 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub 'LU分解 Private Sub decomp(ByRef a, ByRef l, ByRef u) Dim i For i = 0 To N Dim j For j = 0 To N l(i)(j) = 0.0 u(i)(j) = 0.0 Next Next l(0)(0) = 1.0 For j = 0 To N u(0)(j) = a(0)(j) Next For i = 1 To N u(i)(0) = 0.0 l(0)(i) = 0.0 l(i)(0) = a(i)(0) / u(0)(0) Next For i = 1 To N l(i)(i) = 1.0 Dim t: t = a(i)(i) Dim k For k = 0 To i t = t - l(i)(k) * u(k)(i) Next u(i)(i) = t For j = i + 1 To N u(j)(i) = 0.0 l(i)(j) = 0.0 t = a(j)(i) For k = 0 To i t = t - l(j)(k) * u(k)(i) Next l(j)(i) = t / u(i)(i) t = a(i)(j) For k = 0 To i t = t - l(i)(k) * u(k)(j) Next u(i)(j) = t Next Next End Sub '行列の積 Private Sub multiply(ByRef a, ByRef b, ByRef c) Dim i For i = 0 To N Dim j For j = 0 To N Dim s: s = 0.0 Dim k For k = 0 To N s = s + a(i)(k) * b(k)(j) Next c(i)(j) = s Next Next End Sub
Z:\>cscript //nologo Z:\1103.vbs 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
JScript
var N = 4 // LR分解で固有値を求める 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 l = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] var u = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] for (var k = 1; k <= 200; k++) { // LU分解 decomp(a, l, u) // 行列の積 multiply(u, l, a) // 対角要素を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_eigenvalue(a) // 収束判定 var e = 0.0 for (var i = 1; i < N; i++) for (var j = 0; j < i; j++) e += Math.abs(a[i][j]) if (e < 0.00000000001) break } WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") disp_eigenvalue(a) } // LU分解 function decomp(a, l, u) { for (var i = 0; i < N; i++) { for (var j = 0; j < N; j++) { l[i][j] = 0.0 u[i][j] = 0.0 } } l[0][0] = 1.0 for (var j = 0; j < N; j++) u[0][j] = a[0][j] for (var i = 1; i < N; i++) { u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] } for (var i = 1; i < N; i++) { l[i][i] = 1.0 var t = a[i][i] for (var k = 0; k <= i; k++) t -= l[i][k] * u[k][i] u[i][i] = t for (var j = i + 1; j < N; j++) { u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for (var k = 0; k <= i; k++) t -= l[j][k] * u[k][i] l[j][i] = t / u[i][i] t = a[i][j] for (var k = 0; k <= i; k++) t -= l[i][k] * u[k][j] u[i][j] = t } } } // 行列の積 function multiply(a, b, c) { for (var i = 0; i < N; i++) { for (var j = 0; j < N; j++) { var s = 0.0 for (var k = 0; k < N; k++) s += a[i][k] * b[k][j] c[i][j] = s } } } // 対角要素を表示 function disp_eigenvalue(a) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + a[i][i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() }
Z:\>cscript //nologo Z:\1103.js 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 11 9.9993897228 1.0000000001 5.0005095844 2.0001006927 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
PowerShell
set-variable -name N -value 3 -option constant # 対角要素を表示 function disp_eigenvalue($a) { foreach ($i in 0..$N) { Write-Host ([String]::Format("{0,14:F10}`t", $a[$i][$i])) -nonewline } Write-Host } # LU分解 function decomp($a, $l, $u) { foreach ($i in 0..$N) { foreach ($j in 1..$N) { $l[$i][$j] = 0.0 $u[$i][$j] = 0.0 } } $l[0][0] = 1.0 foreach ($j in 0..$N) { $u[0][$j] = $a[0][$j] } foreach ($i in 1..$N) { $u[$i][0] = 0.0 $l[0][$i] = 0.0 $l[$i][0] = $a[$i][0] / $u[0][0] } foreach ($i in 1..$N) { $l[$i][$i] = 1.0 $t = $a[$i][$i] foreach ($k in 0..$i) { $t -= $l[$i][$k] * $u[$k][$i] } $u[$i][$i] = $t if ($i -lt $N) { foreach ($j in ($i + 1)..$N) { $u[$j][$i] = 0 $l[$i][$j] = 0 $t = $a[$j][$i] foreach ($k in 0..$i) { $t -= $l[$j][$k] * $u[$k][$i] } $l[$j][$i] = $t / $u[$i][$i] $t = $a[$i][$j] foreach ($k in 0..$i) { $t -= $l[$i][$k] * $u[$k][$j] } $u[$i][$j] = $t } } } } # 行列の積 function multiply($a, $b, $c) { foreach ($i in 0..$N) { foreach ($j in 0..$N) { $s = 0.0 foreach ($k in 0..$N) { $s += $a[$i][$k] * $b[$k][$j] } $c[$i][$j] = $s } } } # LR分解で固有値を求める $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) $l = (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0) $u = (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0) foreach ($k in 1..200) { # LU分解 decomp $a $l $u # 行列の積 multiply $u $l $a # 対角要素を表示 Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_eigenvalue $a # 収束判定 $e = 0.0 foreach ($i in 1..$N) { foreach ($j in 0..($i - 1)) { $e += [Math]::Abs($a[$i][$j]) } } if ($e -lt 0.00000000001) { break } } Write-Host Write-Host "固有値" disp_eigenvalue $a
Z:\>powershell -file Z:\1103.ps1 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 固有値 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Perl
use constant N => 3; # LR分解で固有値を求める 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 @l, @u; for $k (1..200) { # LU分解 decomp(\@a, \@l, \@u); # 行列の積 multiply(\@u, \@l, \@a); # 対角要素を表示 printf("%3d\t", $k); disp_eigenvalue(\@a); # 収束判定 my $e = 0.0; for $i (1..N) { for $j (0..($i - 1)) { $e += abs($a[$i][$j]); } } last if ($e < 0.00000000001); } print "\neigenvalue\n"; disp_eigenvalue(\@a); } # LU分解 sub decomp { my ($a, $l, $u) = @_; for $i (0..N) { for $j (1..N) { $$l[$i][$j] = 0.0; $$u[$i][$j] = 0.0; } } $$l[0][0] = 1.0; for $j (0..N) { $$u[0][$j] = $$a[0][$j]; } for $i (1..N) { $$u[$i][0] = 0.0; $$l[0][$i] = 0.0; $$l[$i][0] = $$a[$i][0] / $$u[0][0]; } for $i (1..N) { $$l[$i][$i] = 1.0; $t = $$a[$i][$i]; for $k (0..$i) { $t -= $$l[$i][$k] * $$u[$k][$i]; } $$u[$i][$i] = $t; if ($i < N) { for $j (($i + 1)..N) { $$u[$j][$i] = 0; $$l[$i][$j] = 0; $t = $$a[$j][$i]; for $k (0..$i) { $t -= $$l[$j][$k] * $$u[$k][$i]; } $$l[$j][$i] = $t / $$u[$i][$i]; $t = $$a[$i][$j]; for $k (0..$i) { $t -= $$l[$i][$k] * $$u[$k][$j]; } $$u[$i][$j] = $t; } } } } # 行列の積 sub multiply { my ($a, $b, $c) = @_; for $i (0..N) { for $j (0..N) { $s = 0.0; for $k (0..N) { $s += $$a[$i][$k] * $$b[$k][$j]; } $$c[$i][$j] = $s; } } } # 対角要素を表示 sub disp_eigenvalue { my ($a) = @_; for $i (0..N) { printf("%14.10f\t", $$a[$i][$i]); } print("\n"); }
Z:\>perl Z:\1103.pl 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
PHP
<?php define("N", 3); # LR分解で固有値を求める 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]]; $l = [[0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0]]; $u = [[0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0]]; foreach (range(1, 200) as $k) { # LU分解 decomp($a, $l, $u); # 行列の積 multiply($u, $l, $a); # 対角要素を表示 printf("%3d\t", $k); disp_eigenvalue($a); # 収束判定 $e = 0.0; foreach (range(1, N) as $i) { foreach (range(0, ($i - 1)) as $j) { $e += abs($a[$i][$j]); } } if ($e < 0.00000000001) break; } print "\neigenvalue\n"; disp_eigenvalue($a); } # LU分解 function decomp(&$a, &$l, &$u) { foreach (range(0, N) as $i) { foreach (range(1, N) as $j) { $l[$i][$j] = 0.0; $u[$i][$j] = 0.0; } } $l[0][0] = 1.0; foreach (range(0, N) as $j) { $u[0][$j] = $a[0][$j]; } foreach (range(1, N) as $i) { $u[$i][0] = 0.0; $l[0][$i] = 0.0; $l[$i][0] = $a[$i][0] / $u[0][0]; } foreach (range(1, N) as $i) { $l[$i][$i] = 1.0; $t = $a[$i][$i]; foreach (range(0, $i) as $k) { $t -= $l[$i][$k] * $u[$k][$i]; } $u[$i][$i] = $t; if ($i < N) { foreach (range(($i + 1), N) as $j) { $u[$j][$i] = 0; $l[$i][$j] = 0; $t = $a[$j][$i]; foreach (range(0, $i) as $k) { $t -= $l[$j][$k] * $u[$k][$i]; } $l[$j][$i] = $t / $u[$i][$i]; $t = $a[$i][$j]; foreach (range(0, $i) as $k) { $t -= $l[$i][$k] * $u[$k][$j]; } $u[$i][$j] = $t; } } } } # 行列の積 function multiply($a, $b, &$c) { foreach (range(0, N) as $i) { foreach (range(0, N) as $j) { $s = 0.0; foreach (range(0, N) as $k) { $s += $a[$i][$k] * $b[$k][$j]; } $c[$i][$j] = $s; } } } # 対角要素を表示 function disp_eigenvalue($a) { foreach (range(0, N) as $i) { printf("%14.10f\t", $a[$i][$i]); } print("\n"); } ?>
Z:\>php Z:\1103.php 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Python
# coding: Shift_JIS import math N = 4 # 対角要素を表示 def disp_eigenvalue(a): for i in range(0, N, 1): print "%14.10f\t" % a[i][i], print "" # 行列の積 def multiply(a, b, c): for i in range(0, N, 1): for j in range(0, N, 1): s = 0.0 for k in range(0, N, 1): s += a[i][k] * b[k][j] c[i][j] = s # LU分解 def decomp(a, l, u): for i in range(0, N, 1): for j in range(0, N, 1): l[i][j] = 0.0 u[i][j] = 0.0 l[0][0] = 1.0 for j in range(0, N, 1): u[0][j] = a[0][j] for i in range(1, N, 1): u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] for i in range(1, N, 1): l[i][i] = 1.0 t = a[i][i] for k in range(0, i + 1, 1): t -= l[i][k] * u[k][i] u[i][i] = t for j in range(i + 1, N, 1): u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for k in range(0, i + 1, 1): t -= l[j][k] * u[k][i] l[j][i] = t / u[i][i] t = a[i][j] for k in range(0, i + 1, 1): t -= l[i][k] * u[k][j] u[i][j] = t # LR分解で固有値を求める 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]] l = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] u = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] for k in range(1, 200, 1): # LU分解 decomp(a, l, u) # 行列の積 multiply(u, l, a) # 対角要素を表示 print "%3d\t" % k, disp_eigenvalue(a) # 収束判定 e = 0.0 for i in range(1, N, 1): for j in range(0, i, 1): e += abs(a[i][j]) if (e < 0.00000000001): break print "" print "eigenvalue" disp_eigenvalue(a)
Z:\>python Z:\1103.py 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Ruby
# coding: Shift_JIS N = 3 # 対角要素を表示 def disp_eigenvalue(a) (0..N).each do |i| printf("%14.10f\t", a[i][i]) end puts "" end # 行列の積 def multiply(a, b, c) (0..N).each do |i| (0..N).each do |j| s = 0.0 (0..N).each do |k| s += a[i][k] * b[k][j] end c[i][j] = s end end end # LU分解 def decomp(a, l, u) (0..N).each do |i| (0..N).each do |j| l[i][j] = 0.0 u[i][j] = 0.0 end end l[0][0] = 1.0 (0..N).each do |j| u[0][j] = a[0][j] end (1..N).each do |i| u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] end (1..N).each do |i| l[i][i] = 1.0 t = a[i][i] (0..i).each do |k| t -= l[i][k] * u[k][i] end u[i][i] = t ((i + 1)..N).each do |j| u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] (0..i).each do |k| t -= l[j][k] * u[k][i] end l[j][i] = t / u[i][i] t = a[i][j] (0..i).each do |k| t -= l[i][k] * u[k][j] end u[i][j] = t end end end # LR分解で固有値を求める 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]] l = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] u = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] (1..200).each do |k| # LU分解 decomp(a, l, u) # 行列の積 multiply(u, l, a) # 対角要素を表示 printf("%3d\t", k) disp_eigenvalue(a) # 収束判定 e = 0.0 (1..N).each do |i| (0..(i - 1)).each do |j| e += (a[i][j]).abs end end break if (e < 0.00000000001) end puts "" puts "eigenvalue" disp_eigenvalue(a)
Z:\>ruby Z:\1103.rb 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Groovy
N = 3 // LR分解で固有値を求める 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 l = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] as double[][] def u = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] as double[][] for (k in 1..200) { // LU分解 decomp(a, l, u) // 行列の積 multiply(u, l, a) // 対角要素を表示 printf("%3d\t" , k) disp_eigenvalue(a) // 収束判定 e = 0.0 for (i in 1..N) for (j in 0..(i - 1)) e += Math.abs(a[i][j]) if (e < 0.00000000001) break } println() println("eigenvalue") disp_eigenvalue(a) } // LU分解 def decomp(a, l, u) { for (i in 0..N) { for (j in 0..N) { l[i][j] = 0.0 u[i][j] = 0.0 } } l[0][0] = 1.0 for (j in 0..N) u[0][j] = a[0][j] for (i in 1..N) { u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] } for (i in 1..N) { l[i][i] = 1.0 t = a[i][i] for (k in 0..i) t -= l[i][k] * u[k][i] u[i][i] = t if (i < N) { for (j in (i + 1)..N) { u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for (k in 0..i) t -= l[j][k] * u[k][i] l[j][i] = t / u[i][i] t = a[i][j] for (k in 0..i) t -= l[i][k] * u[k][j] u[i][j] = t } } } } // 行列の積 def multiply(a, b, c) { for (i in 0..N) { for (j in 0..N) { s = 0.0 for (k in 0..N) s += a[i][k] * b[k][j] c[i][j] = s } } } // 対角要素を表示 def disp_eigenvalue(a) { for (i in 0..N) printf("%14.10f\t" , a[i][i]) println() }
Z:\>gvy Gvy1103.gvy 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Pascal
program Pas1103(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 対角要素を表示 procedure disp_eigenvalue(a:TwoDimArray); var i:Integer; begin for i := 0 to N do write(format('%14.10f'#9, [a[i][i]])); writeln(); end; // 行列の積 procedure multiply(a:TwoDimArray; b:TwoDimArray; var c:TwoDimArray); var i, j, k:Integer; s:Double; begin for i := 0 to N do begin for j := 0 to N do begin s := 0.0; for k := 0 to N do s := s + (a[i][k] * b[k][j]); c[i][j] := s; end; end; end; // LU分解 procedure decomp(a:TwoDimArray; var l:TwoDimArray; var u:TwoDimArray); var i, j, k:Integer; t:Double; begin for i := 0 to N do begin for j := 0 to N do begin l[i][j] := 0.0; u[i][j] := 0.0; end; end; l[0][0] := 1.0; for j := 0 to N do u[0][j] := a[0][j]; for i := 1 to N do begin u[i][0] := 0.0; l[0][i] := 0.0; l[i][0] := a[i][0] / u[0][0]; end; for i := 1 to N do begin l[i][i] := 1.0; t := a[i][i]; for k := 0 to i do t := t - (l[i][k] * u[k][i]); u[i][i] := t; for j := (i + 1) to N do begin u[j][i] := 0.0; l[i][j] := 0.0; t := a[j][i]; for k := 0 to i do t := t - (l[j][k] * u[k][i]); l[j][i] := t / u[i][i]; t := a[i][j]; for k := 0 to i do t := t - (l[i][k] * u[k][j]); u[i][j] := t; end; end; end; // LR分解で固有値を求める 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)); l :TwoDimArray = ((0.0, 0.0, 0.0, 0.0) ,(0.0, 0.0, 0.0, 0.0) ,(0.0, 0.0, 0.0, 0.0) ,(0.0, 0.0, 0.0, 0.0)); u :TwoDimArray = ((0.0, 0.0, 0.0, 0.0) ,(0.0, 0.0, 0.0, 0.0) ,(0.0, 0.0, 0.0, 0.0) ,(0.0, 0.0, 0.0, 0.0)); e:Double; i, j, k:Integer; begin for k := 1 to 200 do begin // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 write(format('%3d'#9, [k])); disp_eigenvalue(a); // 収束判定 e := 0.0; for i := 1 to N do begin for j := 0 to (i - 1) do e := e + Abs(a[i][j]); end; if e < 0.00000000001 then break; end; writeln(); writeln('eigenvalue'); disp_eigenvalue(a); end.
Z:\>fpc -v0 -l- Pas1103.pp Z:\>Pas1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 141 10.0000000000 5.0000000000 1.9999999982 1.0000000018 142 10.0000000000 5.0000000000 1.9999999991 1.0000000009 143 10.0000000000 5.0000000000 1.9999999995 1.0000000005 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 1.9999999999 1.0000000001 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 150 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.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 Ada1103 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)); l:Long_Float_TwoDimArray; u:Long_Float_TwoDimArray; -- 対角要素を表示 procedure disp_eigenvalue(a:Long_Float_TwoDimArray) is begin for i in 0..N loop Put(a(i, i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; Put_Line(""); end disp_eigenvalue; -- 行列の積 procedure multiply(a:Long_Float_TwoDimArray; b:Long_Float_TwoDimArray; c:in out Long_Float_TwoDimArray) is s:Long_Float; begin for i in 0..N loop for j in 0..N loop s := 0.0; for k in 0..N loop s := s + a(i, k) * b(k, j); end loop; c(i, j) := s; end loop; end loop; end multiply; -- LU分解 procedure decomp(a:Long_Float_TwoDimArray; l:in out Long_Float_TwoDimArray; u:in out Long_Float_TwoDimArray) is t:Long_Float; begin for i in 0..N loop for j in 0..N loop l(i, j) := 0.0; u(i, j) := 0.0; end loop; end loop; l(0, 0) := 1.0; for j in 0..N loop u(0, j) := a(0, j); end loop; for i in 1..N loop u(i, 0) := 0.0; l(0, i) := 0.0; l(i, 0) := a(i, 0) / u(0, 0); end loop; for i in 1..N loop l(i, i) := 1.0; t := a(i, i); for k in 0..i loop t := t - l(i, k) * u(k, i); end loop; u(i, i) := t; for j in (i + 1)..N loop u(j, i) := 0.0; l(i, j) := 0.0; t := a(j, i); for k in 0..i loop t := t - l(j, k) * u(k, i); end loop; l(j, i) := t / u(i, i); t := a(i, j); for k in 0..i loop t := t - l(i, k) * u(k, j); end loop; u(i, j) := t; end loop; end loop; end decomp; -- LR分解で固有値を求める e:Long_Float := 0.0; begin for k in 1..200 loop -- LU分解 decomp(a, l, u); -- 行列の積 multiply(u, l, a); -- 対角要素を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_eigenvalue(a); -- 収束判定 e := 0.0; for i in 1..N loop for j in 0..(i - 1) loop e := e + Abs(a(i, j)); end loop; end loop; if e < 0.00000000001 then exit; end if; end loop; Put_Line(""); Put_Line("eigenvalue"); disp_eigenvalue(a); end Ada1103;
xxxxxx@yyyyyy /Z $ gnatmake Ada1103.adb xxxxxx@yyyyyy /Z $ Ada1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 141 10.0000000000 5.0000000000 1.9999999982 1.0000000018 142 10.0000000000 5.0000000000 1.9999999991 1.0000000009 143 10.0000000000 5.0000000000 1.9999999995 1.0000000005 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 1.9999999999 1.0000000001 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 150 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
VB.NET
Option Explicit Module VB1103 Private Const N As Integer = 3 'LR分解で固有値を求める 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 l(,) As Double = New Double(N,N){} Dim u(,) As Double = New Double(N,N){} For k As Integer= 1 To 200 'LU分解 decomp(a, l, u) '行列の積 multiply(u, l, a) '対角要素を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_eigenvalue(a) '収束判定 Dim e As Double = 0.0 For i As Integer = 1 To N For j As Integer = 0 To i - 1 e += Math.Abs(a(i, j)) Next Next If e < 0.00000000001 Then Exit For Next Console.WriteLine() Console.WriteLine("eigenvalue") disp_eigenvalue(a) End Sub 'LU分解 Private Sub decomp(ByRef a(,) As Double, ByRef l(,) As Double, ByRef u(,) As Double) For i As Integer = 0 To N For j As Integer = 0 To N l(i, j) = 0.0 u(i, j) = 0.0 Next Next l(0, 0) = 1.0 For j As Integer = 0 To N u(0, j) = a(0, j) Next For i As Integer = 1 To N u(i, 0) = 0.0 l(0, i) = 0.0 l(i, 0) = a(i, 0) / u(0, 0) Next For i As Integer = 1 To N l(i, i) = 1.0 Dim t As Double = a(i, i) For k As Integer = 0 To i t -= l(i, k) * u(k, i) Next u(i, i) = t For j As Integer = i + 1 To N u(j, i) = 0.0 l(i, j) = 0.0 t = a(j, i) For k As Integer = 0 To i t -= l(j, k) * u(k, i) Next l(j, i) = t / u(i, i) t = a(i, j) For k As Integer = 0 To i t -= l(i, k) * u(k, j) Next u(i, j) = t Next Next End Sub '行列の積 Private Sub multiply(ByRef a(,) As Double, ByRef b(,) As Double, ByRef c(,) As Double) For i As Integer = 0 To N For j As Integer = 0 To N Dim s As Double = 0.0 For k As Integer = 0 To N s += a(i, k) * b(k, j) Next c(i, j) = s Next Next End Sub '対角要素を表示 Private Sub disp_eigenvalue(ByVal a(,) As Double) For i As Integer = 0 To N Console.Write(String.Format("{0,14:F10}{1}", a(i, i), vbTab)) Next Console.WriteLine() End Sub End Module
Z:\>vbc -nologo VB1103.vb Z:\>VB1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
C#
using System; public class CS1103 { private const int N = 4; // LR分解で固有値を求める 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[,] l = new double[N,N]; double[,] u = new double[N,N]; for (int k = 1; k <= 200; k++) { // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_eigenvalue(a); // 収束判定 double e = 0.0; for (int i = 1; i < N; i++) for (int j = 0; j < i; j++) e += Math.Abs(a[i, j]); if (e < 0.00000000001) break; } Console.WriteLine(); Console.WriteLine("eigenvalue"); disp_eigenvalue(a); } // LU分解 private static void decomp(double[,] a, double[,] l, double[,] u) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { l[i, j] = 0.0; u[i, j] = 0.0; } } l[0, 0] = 1.0; for (int j = 0; j < N; j++) u[0, j] = a[0, j]; for (int i = 1; i < N; i++) { u[i, 0] = 0.0; l[0, i] = 0.0; l[i, 0] = a[i, 0] / u[0, 0]; } for (int i = 1; i < N; i++) { l[i, i] = 1.0; double t = a[i, i]; for (int k = 0; k <= i; k++) t -= l[i, k] * u[k, i]; u[i, i] = t; for (int j = i + 1; j < N; j++) { u[j, i] = 0.0; l[i, j] = 0.0; t = a[j, i]; for (int k = 0; k <= i; k++) t -= l[j, k] * u[k, i]; l[j, i] = t / u[i, i]; t = a[i, j]; for (int k = 0; k <= i; k++) t -= l[i, k] * u[k, j]; u[i, j] = t; } } } // 行列の積 private static void multiply(double[,] a, double[,] b, double[,] c) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double s = 0.0; for (int k = 0; k < N; k++) s += a[i, k] * b[k, j]; c[i,j] = s; } } } // 対角要素を表示 private static void disp_eigenvalue(double[,] a) { for (int i = 0; i < N; i++) Console.Write(string.Format("{0,14:F10}\t", a[i,i])); Console.WriteLine(); } }
Z:\>csc -nologo CS1103.cs Z:\>CS1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Java
import java.lang.*; public class Java1103 { private static final int N = 4; // LR分解で固有値を求める 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[][] l = new double[N][N]; double[][] u = new double[N][N]; for (int k = 1; k <= 200; k++) { // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 System.out.print(String.format("%3d\t", k)); disp_eigenvalue(a); // 収束判定 double e = 0.0; for (int i = 1; i < N; i++) for (int j = 0; j < i; j++) e += Math.abs(a[i][j]); if (e < 0.00000000001) break; } System.out.println(); System.out.println("eigenvalue"); disp_eigenvalue(a); } // LU分解 private static void decomp(double[][] a, double[][] l, double[][] u) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { l[i][j] = 0.0; u[i][j] = 0.0; } } l[0][0] = 1.0; for (int j = 0; j < N; j++) u[0][j] = a[0][j]; for (int i = 1; i < N; i++) { u[i][0] = 0.0; l[0][i] = 0.0; l[i][0] = a[i][0] / u[0][0]; } for (int i = 1; i < N; i++) { l[i][i] = 1.0; double t = a[i][i]; for (int k = 0; k <= i; k++) t -= l[i][k] * u[k][i]; u[i][i] = t; for (int j = i + 1; j < N; j++) { u[j][i] = 0.0; l[i][j] = 0.0; t = a[j][i]; for (int k = 0; k <= i; k++) t -= l[j][k] * u[k][i]; l[j][i] = t / u[i][i]; t = a[i][j]; for (int k = 0; k <= i; k++) t -= l[i][k] * u[k][j]; u[i][j] = t; } } } // 行列の積 private static void multiply(double[][] a, double[][] b, double[][] c) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double s = 0.0; for (int k = 0; k < N; k++) s += a[i][k] * b[k][j]; c[i][j] = s; } } } // 対角要素を表示 private static void disp_eigenvalue(double[][] a) { for (int i = 0; i < N; i++) System.out.print(String.format("%14.10f\t", a[i][i])); System.out.println(); } }
Z:\>javac Java1103.java Z:\>java Java1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
C++
#include <iostream> #include <iomanip> #include <math.h> using namespace std; const int N = 4; // LU分解 void decomp(double a[N][N], double l[N][N], double u[N][N]); // 行列の積 void multiply(double a[N][N], double b[N][N], double c[N][N]); // 対角要素を表示 void disp_eigenvalue(double a[N][N]); // LR分解で固有値を求める 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 l[N][N] = {{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}}; double u[N][N] = {{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}}; for (int k = 1; k <= 200; k++) { // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 cout << setw(3) << k << "\t"; disp_eigenvalue(a); // 収束判定 double e = 0.0; for (int i = 1; i < N; i++) for (int j = 0; j < i; j++) e += fabs(a[i][j]); if (e < 0.00000000001) break; } cout << endl << "eigenvalue" << endl; disp_eigenvalue(a); return 0; } // LU分解 void decomp(double a[N][N], double l[N][N], double u[N][N]) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { l[i][j] = 0.0; u[i][j] = 0.0; } } l[0][0] = 1.0; for (int j = 0; j < N; j++) u[0][j] = a[0][j]; for (int i = 1; i < N; i++) { u[i][0] = 0.0; l[0][i] = 0.0; l[i][0] = a[i][0] / u[0][0]; } for (int i = 1; i < N; i++) { l[i][i] = 1.0; double t = a[i][i]; for (int k = 0; k <= i; k++) t -= l[i][k] * u[k][i]; u[i][i] = t; for (int j = i + 1; j < N; j++) { u[j][i] = 0.0; l[i][j] = 0.0; t = a[j][i]; for (int k = 0; k <= i; k++) t -= l[j][k] * u[k][i]; l[j][i] = t / u[i][i]; t = a[i][j]; for (int k = 0; k <= i; k++) t -= l[i][k] * u[k][j]; u[i][j] = t; } } } // 行列の積 void multiply(double a[N][N], double b[N][N], double c[N][N]) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double s = 0.0; for (int k = 0; k < N; k++) s += a[i][k] * b[k][j]; c[i][j] = s; } } } // 対角要素を表示 void disp_eigenvalue(double a[N][N]) { for (int i = 0; i < N; i++) cout << setw(14) << fixed << setprecision(10) << a[i][i] << "\t"; cout << endl; }
Z:\>bcc32 -q CP1103.cpp cp1103.cpp: Z:\>CP1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 141 10.0000000000 5.0000000000 1.9999999982 1.0000000018 142 10.0000000000 5.0000000000 1.9999999991 1.0000000009 143 10.0000000000 5.0000000000 1.9999999995 1.0000000005 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 1.9999999999 1.0000000001 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 150 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Objective-C
#import <Foundation/Foundation.h> #import <math.h> #define N 4 // LU分解 void decomp(double a[N][N], double l[N][N], double u[N][N]); // 行列の積 void multiply(double a[N][N], double b[N][N], double c[N][N]); // 対角要素を表示 void disp_eigenvalue(double a[N][N]); // LR分解で固有値を求める 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 l[N][N] = {{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}}; double u[N][N] = {{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}}; int k; for (k = 1; k <= 200; k++) { // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 printf("%3d\t", k); disp_eigenvalue(a); // 収束判定 int i, j; double e = 0.0; for (i = 1; i < N; i++) for (j = 0; j < i; j++) e += fabs(a[i][j]); if (e < 0.00000000001) break; } printf("\neigenvalue\n"); disp_eigenvalue(a); return 0; } // LU分解 void decomp(double a[N][N], double l[N][N], double u[N][N]) { int i, j, k; for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { l[i][j] = 0.0; u[i][j] = 0.0; } } l[0][0] = 1.0; for (j = 0; j < N; j++) u[0][j] = a[0][j]; for (i = 1; i < N; i++) { u[i][0] = 0.0; l[0][i] = 0.0; l[i][0] = a[i][0] / u[0][0]; } for (i = 1; i < N; i++) { l[i][i] = 1.0; double t = a[i][i]; for (k = 0; k <= i; k++) t -= l[i][k] * u[k][i]; u[i][i] = t; for (j = i + 1; j < N; j++) { u[j][i] = 0.0; l[i][j] = 0.0; t = a[j][i]; for (k = 0; k <= i; k++) t -= l[j][k] * u[k][i]; l[j][i] = t / u[i][i]; t = a[i][j]; for (k = 0; k <= i; k++) t -= l[i][k] * u[k][j]; u[i][j] = t; } } } // 行列の積 void multiply(double a[N][N], double b[N][N], double c[N][N]) { int i, j, k; for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { double s = 0.0; for (k = 0; k < N; k++) s += a[i][k] * b[k][j]; c[i][j] = s; } } } // 対角要素を表示 void disp_eigenvalue(double a[N][N]) { int i; for (i = 0; i < N; i++) printf("%14.10f\t", a[i][i]); printf("\n"); }
xxxxxx@yyyyyy /Z $ gcc -o OC1103 OC1103.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 141 10.0000000000 5.0000000000 1.9999999982 1.0000000018 142 10.0000000000 5.0000000000 1.9999999991 1.0000000009 143 10.0000000000 5.0000000000 1.9999999995 1.0000000005 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 1.9999999999 1.0000000001 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 150 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
D
import std.stdio; import std.math; const int N = 4; // LR分解で固有値を求める 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 l[N][N] = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]]; double u[N][N] = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]]; foreach (k; 1..200) { // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 writef("%3d\t", k); disp_eigenvalue(a); // 収束判定 double e = 0.0; foreach (i; 1..N) foreach (j; 0..i) e += fabs(a[i][j]); if (e < 0.00000000001) break; } writefln("\neigenvalue"); disp_eigenvalue(a); } // LU分解 void decomp(double[N][N] a, ref double[N][N] l, ref double[N][N] u) { foreach (i; 0..N) { foreach (j; 0..N) { l[i][j] = 0.0; u[i][j] = 0.0; } } l[0][0] = 1.0; foreach (j; 0..N) u[0][j] = a[0][j]; foreach (i; 1..N) { u[i][0] = 0.0; l[0][i] = 0.0; l[i][0] = a[i][0] / u[0][0]; } foreach (i; 1..N) { l[i][i] = 1.0; double t = a[i][i]; foreach (k; 0..i + 1) t -= l[i][k] * u[k][i]; u[i][i] = t; foreach (j; i + 1..N) { u[j][i] = 0.0; l[i][j] = 0.0; t = a[j][i]; foreach (k; 0..i + 1) t -= l[j][k] * u[k][i]; l[j][i] = t / u[i][i]; t = a[i][j]; foreach (k; 0..i + 1) t -= l[i][k] * u[k][j]; u[i][j] = t; } } } // 行列の積 void multiply(double[N][N] a, double[N][N] b, ref double[N][N] c) { foreach (i; 0..N) { foreach (j; 0..N) { double s = 0.0; foreach (k; 0..N) s += a[i][k] * b[k][j]; c[i][j] = s; } } } // 対角要素を表示 void disp_eigenvalue(double[N][N] a) { foreach (i; 0..N) writef("%14.10f\t", a[i][i]); writefln(""); }
Z:\>dmd D1103.d Z:\>D1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 141 10.0000000000 5.0000000000 1.9999999982 1.0000000018 142 10.0000000000 5.0000000000 1.9999999991 1.0000000009 143 10.0000000000 5.0000000000 1.9999999995 1.0000000005 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 1.9999999999 1.0000000001 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 150 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Go
package main import "fmt" import "math" const N = 4 // LR分解で固有値を求める 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 l [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}} var u [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}} for k := 1; k < 200; k++ { // LU分解 decomp(a, &l, &u) // 行列の積 multiply(u, l, &a) // 対角要素を表示 fmt.Printf("%3d\t", k) disp_eigenvalue(a) // 収束判定 var e float64 = 0.0 for i := 1; i < N; i++ { for j := 0; j < i; j++ { e += math.Abs(a[i][j]) } } if (e < 0.00000000001) { break } } fmt.Println("\neigenvector") disp_eigenvalue(a) } // LU分解 func decomp(a[N][N]float64, l *[N][N]float64, u *[N][N]float64) { for i := 0; i < N; i++ { for j := 0; j < N; j++ { l[i][j] = 0.0 u[i][j] = 0.0 } } l[0][0] = 1.0 for j := 0; j < N; j++ { u[0][j] = a[0][j] } for i := 1; i < N; i++ { u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] } for i := 1; i < N; i++ { l[i][i] = 1.0 var t float64 = a[i][i] for k := 0; k <= i; k++ { t -= l[i][k] * u[k][i] } u[i][i] = t for j := i + 1; j < N; j++ { u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for k := 0; k <= i; k++ { t -= l[j][k] * u[k][i] } l[j][i] = t / u[i][i] t = a[i][j] for k := 0; k <= i; k++ { t -= l[i][k] * u[k][j] } u[i][j] = t } } } // 行列の積 func multiply(a[N][N]float64, b[N][N]float64, c *[N][N]float64) { for i := 0; i < N; i++ { for j := 0; j < N; j++ { var s float64 = 0.0 for k := 0; k < N; k++ { s += a[i][k] * b[k][j] } c[i][j] = s } } } // 対角要素を表示 func disp_eigenvalue(a[N][N]float64) { for i := 0; i < N; i++ { fmt.Printf("%14.10f\t", a[i][i]) } fmt.Println("") }
Z:\>go run GO1103.go 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Scala
object Scala1103 { val N = 3 // LR分解で固有値を求める 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 e:Double = 0.0 var k = 1 do { // LU分解 var (l, u) = decomp(a) // 行列の積 a = multiply(u, l) // 対角要素を表示 print("%3d\t".format(k)) disp_eigenvalue(a) k += 1 // 収束判定 e = (for (i <- 0 to N) yield { for (j <- 0 to (i - 1)) yield { Math.abs(a(i)(j)) } }).flatten.sum } while (k <= 200 && e >= 0.00000000001) println() println("eigenvalue") disp_eigenvalue(a) } // LU分解 def decomp(a:Array[Array[Double]]): (Array[Array[Double]], Array[Array[Double]]) = { var l = Array.ofDim[Double](N + 1, N + 1) var u = Array.ofDim[Double](N + 1, N + 1) l(0)(0) = 1.0 for (j <- 0 to N) u(0)(j) = a(0)(j) for (i <- 1 to N) { u(i)(0) = 0.0 l(0)(i) = 0.0 l(i)(0) = a(i)(0) / u(0)(0) } for (i <- 1 to N) { l(i)(i) = 1.0 u(i)(i) = a(i)(i) - (for (k <- 0 to i) yield (l(i)(k) * u(k)(i))).sum for (j <- i + 1 to N) { l(j)(i) = (a(j)(i) - (for (k <- 0 to i) yield (l(j)(k) * u(k)(i))).sum) / u(i)(i) u(i)(j) = a(i)(j) - (for (k <- 0 to i) yield (l(i)(k) * u(k)(j))).sum } } return (l, u) } // 行列の積 def multiply(a:Array[Array[Double]], b:Array[Array[Double]]): Array[Array[Double]] = { (for (i <- 0 to N) yield { (for (j <- 0 to N) yield { (for (k <- 0 to N) yield a(i)(k) * b(k)(j) ).sum }).toArray }).toArray } // 対角要素を表示 def disp_eigenvalue(matrix:Array[Array[Double]]) { (for (i <- 0 to N) yield (matrix(i)(i))).map(col => print("%14.10f\t".format(col))) println() } }
Z:\>scala Scala1103.scala 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 139 10.0000000000 5.0000000000 2.0000000032 0.9999999968 140 10.0000000000 5.0000000000 2.0000000016 0.9999999984 141 10.0000000000 5.0000000000 2.0000000008 0.9999999992 142 10.0000000000 5.0000000000 2.0000000004 0.9999999996 143 10.0000000000 5.0000000000 2.0000000002 0.9999999998 144 10.0000000000 5.0000000000 2.0000000001 0.9999999999 145 10.0000000000 5.0000000000 2.0000000000 1.0000000000 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
F#
module Fs1103 open System let N = 3 // LU分解 let decomp (a:float[][]) : float[][] * float[][] = let l:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] let u:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] l.[0].[0] <- 1.0 for j in [0..N] do u.[0].[j] <- a.[0].[j] for i in [1..N] do u.[i].[0] <- 0.0 l.[0].[i] <- 0.0 l.[i].[0] <- a.[i].[0] / u.[0].[0] for i in [1..N] do l.[i].[i] <- 1.0 u.[i].[i] <- a.[i].[i] - ([| for k in 0..i -> l.[i].[k] * u.[k].[i] |] |> Array.sum) for j in [(i + 1)..N] do l.[j].[i] <- (a.[j].[i] - ([| for k in 0..i -> l.[j].[k] * u.[k].[i] |] |> Array.sum)) / u.[i].[i] u.[i].[j] <- a.[i].[j] - ([| for k in 0..i -> l.[i].[k] * u.[k].[j] |] |> Array.sum) (l, u) // 行列の積 let multiply (a:float[][]) (b:float[][]) : float[][] = [| for i in 0..N -> [| for j in 0..N -> [| for k in 0..N -> a.[i].[k] * b.[k].[j] |] |> Array.sum |] |] // 対角要素を表示 let disp_eigenvalue (a:float[][]) = for i in [0..N] do printf "%14.10f\t" a.[i].[i] printfn "" // LR分解で固有値を求める 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 e = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 200 && e >= 0.00000000001) do if sw then sw <- false // LU分解 let (l, u) = decomp a // 行列の積 a <- multiply u l // 対角要素を表示 printf "%3d\t" k disp_eigenvalue a k <- k + 1 // 収束判定 e <- [| for i in 0..N -> [| for j in 0..(i - 1) -> Math.Abs(a.[i].[j]) |] |] |> Array.collect id |> Array.sum printfn "" printfn "eigenvalue" disp_eigenvalue a exit 0
Z:\>fsi --nologo --quiet Fs1103.fs 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 省略 139 10.0000000000 5.0000000000 2.0000000032 0.9999999968 140 10.0000000000 5.0000000000 2.0000000016 0.9999999984 141 10.0000000000 5.0000000000 2.0000000008 0.9999999992 142 10.0000000000 5.0000000000 2.0000000004 0.9999999996 143 10.0000000000 5.0000000000 2.0000000002 0.9999999998 144 10.0000000000 5.0000000000 2.0000000001 0.9999999999 145 10.0000000000 5.0000000000 2.0000000000 1.0000000000 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000