さまざまな言語で数値計算
Only Do What Only You Can Do
QR分解法
QR分解を使って, 全ての固有値を求める
対称行列 $ \boldsymbol{A} $ を 正規直交行列 $ \boldsymbol{Q} $ と 右三角行列 $ \boldsymbol{R} $ の積に分解し,
- $ \boldsymbol{A}_k = \boldsymbol{Q}_k \boldsymbol{R}_k $
- $ \boldsymbol{A}_{k+1} = \boldsymbol{R}_k \boldsymbol{Q}_k $
VBScript
Option Explicit Private Const N = 3 'QR分解で固有値を求める 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 q: q = 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 r: r = 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 'QR分解 Call decomp(a, q, r) '行列の積 Call multiply(r, q, 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 'QR分解 Private Sub decomp(ByRef a, ByRef q, ByRef r) Dim x: x = Array(0.0, 0.0, 0.0, 0.0) Dim k For k = 0 To N Dim i For i = 0 To N x(i) = a(i)( k) Next Dim j For j = 0 To k - 1 Dim t: t = 0 For i = 0 To N t = t + a(i)(k) * q(i)(j) Next r(j)(k) = t r(k)(j) = 0 For i = 0 To N x(i) = x(i) - t * q(i)(j) Next Next t = 0 For i = 0 To N t = t + x(i) * x(i) Next r(k)(k) = Sqr(t) For i = 0 To N q(i)(k) = x(i) / r(k)(k) 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:\1104.vbs 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
JScript
var N = 4 // QR分解で固有値を求める 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 q = [[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 r = [[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++) { // QR分解 decomp(a, q, r) // 行列の積 multiply(r, q, 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) } // QR分解 function decomp(a, q, r) { var x = [0.0 ,0.0 ,0.0 ,0.0] for (var k = 0; k < N; k++) { for (var i = 0; i < N; i++) x[i] = a[i][k] for (var j = 0; j < k; j++) { var t = 0.0 for (var i = 0; i < N; i++) t += a[i][k] * q[i][j] r[j][k] = t r[k][j] = 0.0 for (var i = 0; i < N; i++) x[i] -= t * q[i][j] } var s = 0.0 for (var i = 0; i < N; i++) s += x[i] * x[i] r[k][k] = Math.sqrt(s) for (var i = 0; i < N; i++) q[i][k] = x[i] / r[k][k] } } // 行列の積 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:\1104.js 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 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 } # QR分解 function decomp($a, $q, $r) { $x = (0.0, 0.0, 0.0, 0.0) foreach ($k in 0..$N) { foreach ($i in 0..$N) { $x[$i] = $a[$i][$k] } if ($k -gt 0) { foreach ($j in 0..($k - 1)) { $t = 0 foreach ($i in 0..$N) { $t += $a[$i][$k] * $q[$i][$j] } $r[$j][$k] = $t $r[$k][$j] = 0 foreach ($i in 0..$N) { $x[$i] -= $t * $q[$i][$j] } } } $t = 0 foreach ($i in 0..$N) { $t += $x[$i] * $x[$i] } $r[$k][$k] = [Math]::Sqrt($t) foreach ($i in 0..$N) { $q[$i][$k] = $x[$i] / $r[$k][$k] } } } # 行列の積 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 } } } # QR分解で固有値を求める $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) $q = (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) $r = (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) { # QR分解 decomp $a $q $r # 行列の積 multiply $r $q $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:\1104.ps1 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 固有値 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Perl
use constant N => 3; # QR分解で固有値を求める 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 @q, @r; for $k (1..100) { # QR分解 decomp(\@a, \@q, \@r); # 行列の積 multiply(\@r, \@q, \@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); } # QR分解 sub decomp { my ($a, $q, $r) = @_; my @x; for $k (0..N) { for $i (0..N) { $x[$i] = $$a[$i][$k]; } for $j (0..($k - 1)) { my $t = 0.0; for $i (0..N) { $t += $$a[$i][$k] * $$q[$i][$j]; } $$r[$j][$k] = $t; $$r[$k][$j] = 0.0; for $i (0..N) { $x[$i] -= $t * $$q[$i][$j]; } } my $s = 0.0; for $i (0..N) { $s += $x[$i] * $x[$i]; } $$r[$k][$k] = sqrt($s); for $i (0..N) { $$q[$i][$k] = $x[$i] / $$r[$k][$k]; } } } # 行列の積 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:\1104.pl 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
PHP
<?php define("N", 3); # QR分解で固有値を求める 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]]; $q = [[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]]; $r = [[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) { # QR分解 decomp($a, $q, $r); # 行列の積 multiply($r, $q, $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); } # QR分解 function decomp($a, &$q, &$r) { $x = [0.0, 0.0, 0.0, 0.0]; foreach (range(0, N) as $k) { foreach (range(0, N) as $i) { $x[$i] = $a[$i][$k]; } if ($k > 0) { foreach (range(0, ($k - 1)) as $j) { $t = 0.0; foreach (range(0, N) as $i) { $t += $a[$i][$k] * $q[$i][$j]; } $r[$j][$k] = $t; $r[$k][$j] = 0.0; foreach (range(0, N) as $i) { $x[$i] -= $t * $q[$i][$j]; } } } $s = 0.0; foreach (range(0, N) as $i) { $s += $x[$i] * $x[$i]; } $r[$k][$k] = sqrt($s); foreach (range(0, N) as $i) { $q[$i][$k] = $x[$i] / $r[$k][$k]; } } } # 行列の積 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:\1104.php 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 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 # QR分解 def decomp(a, q, r): x = [0.0 ,0.0 ,0.0 ,0.0] for k in range(0, N, 1): for i in range(0, N, 1): x[i] = a[i][k] for j in range(0, k, 1): t = 0.0 for i in range(0, N, 1): t += a[i][k] * q[i][j] r[j][k] = t r[k][j] = 0.0 for i in range(0, N, 1): x[i] -= t * q[i][j] s = 0.0 for i in range(0, N, 1): s += x[i] * x[i] r[k][k] = math.sqrt(s) for i in range(0, N, 1): q[i][k] = x[i] / r[k][k] # QR分解で固有値を求める 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]] q = [[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]] r = [[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): # QR分解 decomp(a, q, r) # 行列の積 multiply(r, q, 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:\1104.py 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 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 # QR分解 def decomp(a, q, r) x = [0.0 ,0.0 ,0.0 ,0.0] (0..N).each do |k| (0..N).each do |i| x[i] = a[i][k] end (0..(k - 1)).each do |j| t = 0.0 (0..N).each do |i| t += a[i][k] * q[i][j] end r[j][k] = t r[k][j] = 0.0 (0..N).each do |i| x[i] -= t * q[i][j] end end s = 0.0 (0..N).each do |i| s += x[i] * x[i] end r[k][k] = Math.sqrt(s) (0..N).each do |i| q[i][k] = x[i] / r[k][k] end end end # QR分解で固有値を求める 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]] q = [[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]] r = [[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| # QR分解 decomp(a, q, r) # 行列の積 multiply(r, q, 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:\1104.rb 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Groovy
N = 3 // QR分解で固有値を求める 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 q = [[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 r = [[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) { // QR分解 decomp(a, q, r) // 行列の積 multiply(r, q, 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) } // QR分解 def decomp(a, q, r) { x = [0.0 ,0.0 ,0.0 ,0.0] as double[] for (k in 0..N) { for (i in 0..N) x[i] = a[i][k] if (k > 0) { for (j in 0..(k - 1)) { t = 0.0 for (i in 0..N) t += a[i][k] * q[i][j] r[j][k] = t r[k][j] = 0.0 for (i in 0..N) x[i] -= t * q[i][j] } } s = 0.0 for (i in 0..N) s += x[i] * x[i] r[k][k] = Math.sqrt(s) for (i in 0..N) q[i][k] = x[i] / r[k][k] } } // 行列の積 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 Gvy1104.gvy 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Pascal
program Pas1104(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; // QR分解 procedure decomp(a:TwoDimArray; var q:TwoDimArray; var r:TwoDimArray); var i, j, k:Integer; t, s:Double; x :array [0..N] of Double; begin for k := 0 to N do begin for i := 0 to N do x[i] := a[i][k]; for j := 0 to (k - 1) do begin t := 0.0; for i := 0 to N do t := t + a[i][k] * q[i][j]; r[j][k] := t; r[k][j] := 0.0; for i := 0 to N do x[i] := x[i] - t * q[i][j]; end; s := 0.0; for i := 0 to N do s := s + x[i] * x[i]; r[k][k] := Sqrt(s); for i := 0 to N do q[i][k] := x[i] / r[k][k]; end; end; // QR分解で固有値を求める 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)); q :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)); r :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 // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, 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- Pas1104.pp Z:\>Pas1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 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 Ada1104 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)); q:Long_Float_TwoDimArray; r: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; -- QR分解 procedure decomp(a:Long_Float_TwoDimArray; q:in out Long_Float_TwoDimArray; r:in out Long_Float_TwoDimArray) is t, s:Long_Float; x:Long_Float_Array; begin for k in 0..N loop for i in 0..N loop x(i) := a(i, k); end loop; for j in 0..(k - 1) loop t := 0.0; for i in 0..N loop t := t + a(i, k) * q(i, j); end loop; r(j, k) := t; r(k, j) := 0.0; for i in 0..N loop x(i) := x(i) - t * q(i, j); end loop; end loop; s := 0.0; for i in 0..N loop s := s + x(i) * x(i); end loop; r(k, k) := Sqrt(s); for i in 0..N loop q(i, k) := x(i) / r(k, k); end loop; end loop; end decomp; -- QR分解で固有値を求める e:Long_Float := 0.0; begin for k in 1..200 loop -- QR分解 decomp(a, q, r); -- 行列の積 multiply(r, q, 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 Ada1104;
xxxxxx@yyyyyy /Z $ gnatmake Ada1104.adb xxxxxx@yyyyyy /Z $ Ada1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
VB.NET
Option Explicit Module VB1104 Private Const N As Integer = 3 'QR分解で固有値を求める 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 q(,) As Double = New Double(N,N){} Dim r(,) As Double = New Double(N,N){} For k As Integer= 1 To 100 'QR分解 decomp(a, q, r) '行列の積 multiply(r, q, 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 'QR分解 Private Sub decomp(ByRef a(,) As Double, ByRef q(,) As Double, ByRef r(,) As Double) Dim x() As Double = New Double(N){} For k As Integer = 0 To N For i As Integer = 0 To N x(i) = a(i, k) Next For j As Integer = 0 To k - 1 Dim t As Double = 0.0 For i As Integer = 0 To N t += a(i, k) * q(i, j) Next r(j, k) = t r(k, j) = 0.0 For i As Integer = 0 To N x(i) -= t * q(i, j) Next Next Dim s As Double = 0.0 For i As Integer = 0 To N s += x(i) * x(i) Next r(k, k) = Math.Sqrt(s) For i As Integer = 0 To N q(i, k) = x(i) / r(k, k) 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 VB1104.vb Z:\>VB1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
C#
using System; public class CS1104 { private const int N = 4; // QR分解で固有値を求める 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[,] q = new double[N,N]; double[,] r = new double[N,N]; for (int k = 1; k <= 100; k++) { // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, a); // 対角要素を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_eigenvalue(a); // 収束判定 double e = 0.0; for (int i = 0; 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); } // QR分解 private static void decomp(double[,] a, double[,] q, double[,] r) { double[] x = new double[N]; for (int k = 0; k < N; k++) { for (int i = 0; i < N; i++) x[i] = a[i, k]; for (int j = 0; j < k; j++) { double t = 0.0; for (int i = 0; i < N; i++) t += a[i, k] * q[i, j]; r[j, k] = t; r[k, j] = 0.0; for (int i = 0; i < N; i++) x[i] -= t * q[i, j]; } double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; r[k,k] = Math.Sqrt(s); for (int i = 0; i < N; i++) q[i, k] = x[i] / r[k, k]; } } // 行列の積 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 CS1104.cs Z:\>CS1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Java
import java.lang.*; public class Java1104 { private static final int N = 4; // QR分解で固有値を求める 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[][] q = new double[N][N]; double[][] r = new double[N][N]; for (int k = 1; k <= 200; k++) { // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, 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("固有値"); disp_eigenvalue(a); } // QR分解 private static void decomp(double[][] a, double[][] q, double[][] r) { double[] x = new double[N]; for (int k = 0; k < N; k++) { for (int i = 0; i < N; i++) x[i] = a[i][k]; for (int j = 0; j < k; j++) { double t = 0.0; for (int i = 0; i < N; i++) t += a[i][k] * q[i][j]; r[j][k] = t; r[k][j] = 0.0; for (int i = 0; i < N; i++) x[i] -= t * q[i][j]; } double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; r[k][k] = Math.sqrt(s); for (int i = 0; i < N; i++) q[i][k] = x[i] / r[k][k]; } } // 行列の積 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 Java1104.java Z:\>java Java1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 固有値 10.0000000000 5.0000000000 2.0000000000 1.0000000000
C++
#include <iostream> #include <iomanip> #include <math.h> using namespace std; const int N = 4; // QR分解 void decomp(double a[N][N], double q[N][N], double r[N][N]); // 行列の積 void multiply(double a[N][N], double b[N][N], double c[N][N]); // 対角要素を表示 void disp_eigenvalue(double a[N][N]); // QR分解で固有値を求める 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 q[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 r[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++) { // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, 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; } // QR分解 void decomp(double a[N][N], double q[N][N], double r[N][N]) { double x[N] = {0.0 ,0.0 ,0.0 ,0.0}; for (int k = 0; k < N; k++) { for (int i = 0; i < N; i++) x[i] = a[i][k]; for (int j = 0; j < k; j++) { double t = 0.0; for (int i = 0; i < N; i++) t += a[i][k] * q[i][j]; r[j][k] = t; r[k][j] = 0.0; for (int i = 0; i < N; i++) x[i] -= t * q[i][j]; } double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; r[k][k] = sqrt(s); for (int i = 0; i < N; i++) q[i][k] = x[i] / r[k][k]; } } // 行列の積 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 CP1104.cpp cp1104.cpp: Z:\>CP1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 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 // QR分解 void decomp(double a[N][N], double q[N][N], double r[N][N]); // 行列の積 void multiply(double a[N][N], double b[N][N], double c[N][N]); // 対角要素を表示 void disp_eigenvalue(double a[N][N]); // QR分解で固有値を求める 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 q[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 r[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++) { // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, a); // 対角要素を表示 printf("%3d\t", k); disp_eigenvalue(a); // 収束判定 int i, j; double e = 0.0; for (i = 0; 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; } // QR分解 void decomp(double a[N][N], double q[N][N], double r[N][N]) { int i, j, k; double x[N] = {0.0 ,0.0 ,0.0 ,0.0}; for (k = 0; k < N; k++) { for (i = 0; i < N; i++) x[i] = a[i][k]; for (j = 0; j < k; j++) { double t = 0.0; for (i = 0; i < N; i++) t += a[i][k] * q[i][j]; r[j][k] = t; r[k][j] = 0.0; for (i = 0; i < N; i++) x[i] -= t * q[i][j]; } double s = 0.0; for (i = 0; i < N; i++) s += x[i] * x[i]; r[k][k] = sqrt(s); for (i = 0; i < N; i++) q[i][k] = x[i] / r[k][k]; } } // 行列の積 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 OC1104 OC1104.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 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; // QR分解で固有値を求める 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 q[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 r[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) { // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, 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); } // QR分解 void decomp(double[N][N] a, ref double[N][N] q, ref double[N][N] r) { double x[N] = [0.0 ,0.0 ,0.0 ,0.0]; foreach (k; 0..N) { foreach (i; 0..N) x[i] = a[i][k]; foreach (j; 0..k) { double t = 0.0; foreach (i; 0..N) t += a[i][k] * q[i][j]; r[j][k] = t; r[k][j] = 0.0; foreach (i; 0..N) x[i] -= t * q[i][j]; } double s = 0.0; foreach (i; 0..N) s += x[i] * x[i]; r[k][k] = sqrt(s); foreach (i; 0..N) q[i][k] = x[i] / r[k][k]; } } // 行列の積 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 D1104.d Z:\>D1104 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 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 // QR分解で固有値を求める 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 q [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 r [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++ { // QR分解 decomp(a, &q, &r) // 行列の積 multiply(r, q, &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("\neigenvalue") disp_eigenvalue(a) } // QR分解 func decomp(a[N][N]float64, q *[N][N]float64, r *[N][N]float64) { var x = []float64 {0.0 ,0.0 ,0.0 ,0.0} for k := 0; k < N; k++ { for i := 0; i < N; i++ { x[i] = a[i][k] } for j := 0; j < k; j++ { var t float64 = 0.0 for i := 0; i < N; i++ { t += a[i][k] * q[i][j] } r[j][k] = t r[k][j] = 0.0 for i := 0; i < N; i++ { x[i] -= t * q[i][j] } } var s float64 = 0.0 for i := 0; i < N; i++ { s += x[i] * x[i] } r[k][k] = math.Sqrt(s) for i := 0; i < N; i++ { q[i][k] = x[i] / r[k][k] } } } // 行列の積 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 GO1104.go 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
Scala
object Scala1104 { val N = 3 // QR分解で固有値を求める 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 { // QR分解 var (q, r) = decomp(a) // 行列の積 a = multiply(r, q) // 対角要素を表示 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) } // QR分解 def decomp(a:Array[Array[Double]]): (Array[Array[Double]], Array[Array[Double]]) = { var x = new Array[Double](N + 1) var q = Array.ofDim[Double](N + 1, N + 1) var r = Array.ofDim[Double](N + 1, N + 1) for (k <- 0 to N) { for (i <- 0 to N) x(i) = a(i)(k) for (j <- 0 to k - 1) { r(j)(k) = (for (i <- 0 to N) yield (a(i)(k) * q(i)(j))).sum r(k)(j) = 0.0 for (i <- 0 to N) x(i) -= r(j)(k) * q(i)(j) } r(k)(k) = Math.sqrt(x.map(n => n * n).sum) for (i <- 0 to N) q(i)(k) = x(i) / r(k)(k) } return (q, r) } // 行列の積 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 Scala1104.scala 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000
F#
module Fs1104 open System let N = 3 // QR分解 let decomp (a:float[][]) : float[][] * float[][] = let q:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] let r:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] for k in [0..N] do let x:float[] = [| for i in 0..N -> a.[i].[k] |] for j in [0..(k - 1)] do r.[j].[k] <- [| for i in 0..N -> a.[i].[k] * q.[i].[j] |] |> Array.sum r.[k].[j] <- 0.0 for i in [0..N] do x.[i] <- x.[i] - r.[j].[k] * q.[i].[j] r.[k].[k] <- Math.Sqrt(x |> Array.map(fun n -> n * n) |> Array.sum) for i in [0..N] do q.[i].[k] <- x.[i] / r.[k].[k] (q, r) // 行列の積 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 "" // QR分解で固有値を求める 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 // QR分解 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 Fs1104.fs 1 9.6046511628 1.1012311902 4.8997514499 2.3943661972 2 9.9219788334 1.0010980897 5.0142271522 2.0626959248 3 9.9805335651 1.0000111820 5.0095550207 2.0099002322 4 9.9951218389 1.0000001123 5.0033019159 2.0015761328 5 9.9987795937 1.0000000011 5.0009686042 2.0002518010 6 9.9996948428 1.0000000000 5.0002648858 2.0000402713 7 9.9999237072 1.0000000000 5.0000698501 2.0000064427 8 9.9999809266 1.0000000000 5.0000180426 2.0000010308 9 9.9999952316 1.0000000000 5.0000046034 2.0000001649 10 9.9999988079 1.0000000000 5.0000011657 2.0000000264 省略 85 10.0000000000 5.0000000000 2.0000000000 1.0000000000 86 10.0000000000 5.0000000000 2.0000000000 1.0000000000 87 10.0000000000 5.0000000000 2.0000000000 1.0000000000 88 10.0000000000 5.0000000000 2.0000000000 1.0000000000 89 10.0000000000 5.0000000000 2.0000000000 1.0000000000 90 10.0000000000 5.0000000000 2.0000000000 1.0000000000 91 10.0000000000 5.0000000000 2.0000000000 1.0000000000 92 10.0000000000 5.0000000000 2.0000000000 1.0000000000 93 10.0000000000 5.0000000000 2.0000000000 1.0000000000 94 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000