さまざまな言語で数値計算
Only Do What Only You Can Do
ヤコビ法
ヤコビ法を使って, 全ての固有値と, 対応する固有ベクトルを求める
対称行列 $ \boldsymbol{A} $ の非対角成分のうち絶対値が最大の成分を $ a_{ij} $ とするとき,
$ \boldsymbol{A} $ に $ \boldsymbol{P} $ と, 転置行列 $ \boldsymbol{P}^T $ を左右からかけてできた行列 $ \boldsymbol{B} = \boldsymbol{P}^T\boldsymbol{AP} $ の固有値と固有ベクトルは 元の行列 $ \boldsymbol{A} $ から変化しない. これを相似変換という.
$ b_{ij}=b_{ji}=0 $ となるように $ \theta $ の値を設定して処理を反復すると, 固有値が対角成分に並んだ対角行列に収束する.
VBScript
Option Explicit Private Const N = 3 Private Const PI = 3.14159265358979 'ヤコビ法で固有値を求める 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 v: v = Array(_ Array(1.0 ,0.0 ,0.0 ,0.0), _ Array(0.0 ,1.0 ,0.0 ,0.0), _ Array(0.0 ,0.0 ,1.0 ,0.0), _ Array(0.0 ,0.0 ,0.0 ,1.0)) 'ヤコビ法 Call jacobi(a, v) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" Call disp_eigenvalue(a) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvector" Call disp_eigenvector(v) End Sub 'ヤコビ法 Private Sub jacobi(ByRef a, ByRef v) Dim i, j, k Dim max_val Dim p, q For k = 1 To 100 '最大値を探す max_val = 0.0 For i = 0 To N For j = i + 1 To N If (max_val < Abs(a(i)(j))) Then max_val = Abs(a(i)(j)) p = i q = j End If Next Next 'θ を求める Dim t If Abs(a(p)(p) - a(q)(q)) < 0.00000000001 Then 'a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = PI / 4.0 If (a(p)(p) < 0) Then t = -t End If Else 'a_{pp} ≠ a_{qq} のとき t = Atn(2.0 * a(p)(q) / (a(p)(p) - a(q)(q))) / 2.0 End If 'θ を使って 行列 U を作成し、A = U^t × A × U Dim c: c = Cos(t) Dim s: s = Sin(t) 'U^t × A Dim t1, t2 For i = 0 To N t1 = a(p)(i) * c + a(q)(i) * s t2 = -a(p)(i) * s + a(q)(i) * c a(p)(i) = t1 a(q)(i) = t2 '固有ベクトル t1 = v(p)(i) * c + v(q)(i) * s t2 = -v(p)(i) * s + v(q)(i) * c v(p)(i) = t1 v(q)(i) = t2 Next 'A × U For i = 0 To N t1 = a(i)(p) * c + a(i)(q) * s t2 = -a(i)(p) * s + a(i)(q) * c a(i)(p) = t1 a(i)(q) = t2 Next '行列の対角要素を表示 (固有値) WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_eigenvalue(a) '収束判定 If max_val < 0.00000000001 Then Exit For Next End Sub '行列の対角要素を表示 Private Sub disp_eigenvalue(ByVal x) Dim i For i = 0 To N WScript.StdOut.Write Right(Space(14) & FormatNumber(x(i)(i), 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '固有ベクトルを表示 Private Sub disp_eigenvector(ByVal matrix) Dim row, col For Each row In matrix normarize(row) disp_vector(row) WScript.StdOut.WriteLine Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x()) Dim s: s = 0.0 Dim i For i = 0 To N s = s + x(i) * x(i) Next s = Sqr(s) For i = 0 To N x(i) = x(i) / s Next End Sub
Z:\>cscript //nologo Z:\1105.vbs 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
JScript
var N = 4 // ヤコビ法で固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var v = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]] // ヤコビ法 jacobi(a, v) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") disp_eigenvalue(a) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvector") disp_eigenvector(v) } // ヤコビ法 function jacobi(a, v) { for (var k = 1; k <= 100; k++) { // 最大値を探す var p = 0 var q = 0 var max_val = 0.0 for (var i = 0; i < N; i++) { for (var j = i + 1; j < N; j++) { if (max_val < Math.abs(a[i][j])) { max_val = Math.abs(a[i][j]) p = i q = j } } } // θ を求める var t = 0.0 if (Math.abs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 if (a[p][p] < 0) t = -t } else { // a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 } // θ を使って 行列 U を作成し、A = U^t × A × U var c = Math.cos(t) var s = Math.sin(t) // U^t × A var t1 = 0.0 var t2 = 0.0 for (var i = 0; i < N; i++) { t1 = a[p][i] * c + a[q][i] * s t2 = -a[p][i] * s + a[q][i] * c a[p][i] = t1 a[q][i] = t2 // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s t2 = -v[p][i] * s + v[q][i] * c v[p][i] = t1 v[q][i] = t2 } // A × U for (var i = 0; i < N; i++) { t1 = a[i][p] * c + a[i][q] * s t2 = -a[i][p] * s + a[i][q] * c a[i][p] = t1 a[i][q] = t2 } // 対角要素を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_eigenvalue(a) // 収束判定 if (max_val < 0.00000000001) break } } // 対角要素を表示 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() } // 固有ベクトルを表示 function disp_eigenvector(matrix) { for (var i = 0; i < N; i++) { var row = [0.0 ,0.0 ,0.0 ,0.0] for (var j = 0; j < N; j++) row[j] = matrix[i][j] normarize(row) disp_vector(row) } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 正規化 (ベクトルの長さを1にする) function normarize(x) { var s = 0.0 for (var i = 0; i < N; i++) s += x[i] * x[i] s = Math.sqrt(s) for (var i = 0; i < N; i++) x[i] /= s }
Z:\>cscript //nologo Z:\1105.js 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
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 } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 正規化 (ベクトルの長さを1にする) function normarize($x) { $s = 0.0 foreach ($i in 0..$N) { $s += $x[$i] * $x[$i] } $s = [Math]::Sqrt($s) foreach ($i in 0..$N) { $x[$i] /= $s } } # 固有ベクトルを表示 function disp_eigenvector($matrix) { foreach ($row in $matrix) { normarize $row disp_vector $row } } # ヤコビ法 function jacobi($a, $v) { foreach ($k in 1..100) { # 最大値を探す $max_val = 0.0 foreach ($i in 0..($N - 1)) { foreach ($j in ($i + 1)..$N) { if ($max_val -lt [Math]::Abs($a[$i][$j])) { $max_val = [Math]::Abs($a[$i][$j]) $p = $i $q = $j } } } # θ を求める if ([Math]::Abs($a[$p][$p] - $a[$q][$q]) -lt 0.00000000001) { # a_{pp} = a_{qq} のとき、回転角tをπ/4にする $t = [Math]::PI / 4.0 if ($a[$p][$p] -lt 0.0) { $t = -$t } } else { # a_{pp} ≠ a_{qq} のとき $t = [Math]::Atan(2.0 * $a[$p][$q] / ($a[$p][$p] - $a[$q][$q])) / 2.0 } # θ を使って 行列 U を作成し、A = U^t × A × U $c = [Math]::Cos($t) $s = [Math]::Sin($t) # U^t × A foreach ($i in 0..$N) { $t1 = $a[$p][$i] * $c + $a[$q][$i] * $s $t2 = -$a[$p][$i] * $s + $a[$q][$i] * $c $a[$p][$i] = $t1 $a[$q][$i] = $t2 # 固有ベクトル $t1 = $v[$p][$i] * $c + $v[$q][$i] * $s $t2 = -$v[$p][$i] * $s + $v[$q][$i] * $c $v[$p][$i] = $t1 $v[$q][$i] = $t2 } # A × U foreach ($i in 0..$N) { $t1 = $a[$i][$p] * $c + $a[$i][$q] * $s $t2 = -$a[$i][$p] * $s + $a[$i][$q] * $c $a[$i][$p] = $t1 $a[$i][$q] = $t2 } # 行列の対角要素を表示 (固有値) Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_eigenvalue $a # 収束判定 if ($max_val -lt 0.00000000001) { break } } } # ヤコビ法で固有値を求める $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) $v = (1.0 ,0.0 ,0.0 ,0.0), (0.0 ,1.0 ,0.0 ,0.0), (0.0 ,0.0 ,1.0 ,0.0), (0.0 ,0.0 ,0.0 ,1.0) # ヤコビ法 jacobi $a $v Write-Host Write-Host "固有値" disp_eigenvalue $a Write-Host Write-Host "固有ベクトル" disp_eigenvector $v
Z:\>powershell -file Z:\1105.ps1 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 固有値 10.0000000000 1.0000000000 5.0000000000 2.0000000000 固有ベクトル 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
Perl
use Math::Trig; use constant N => 3; # ヤコビ法で固有値を求める main(); sub main { my @a = ([5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]); my @v = ([1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]); # ヤコビ法 jacobi(\@a, \@v); print "\neigenvalue\n"; disp_eigenvalue(\@a); print "\neigenvector\n"; disp_eigenvector(\@v); } # ヤコビ法 sub jacobi { my ($a, $v) = @_; for $k (1..100) { # 最大値を探す $max_val = 0.0; for $i (0..(N - 1)) { for $j (($i + 1)..N) { if ($max_val < abs($$a[$i][$j])) { $max_val = abs($$a[$i][$j]); $p = $i; $q = $j; } } } # θ を求める if (abs($$a[$p][$p] - $$a[$q][$q]) < 0.00000000001) { # a_{pp} = a_{qq} のとき、回転角tをπ/4にする $t = pi / 4.0; $t = -$t if ($$a[$p][$p] < 0.0) } else { # a_{pp} ≠ a_{qq} のとき $t = atan(2.0 * $$a[$p][$q] / ($$a[$p][$p] - $$a[$q][$q])) / 2.0; } # θ を使って 行列 U を作成し、A = U^t × A × U $c = cos($t); $s = sin($t); # U^t × A for $i (0..N) { $t1 = $$a[$p][$i] * $c + $$a[$q][$i] * $s; $t2 = -$$a[$p][$i] * $s + $$a[$q][$i] * $c; $$a[$p][$i] = $t1; $$a[$q][$i] = $t2; # 固有ベクトル $t1 = $$v[$p][$i] * $c + $$v[$q][$i] * $s; $t2 = -$$v[$p][$i] * $s + $$v[$q][$i] * $c; $$v[$p][$i] = $t1; $$v[$q][$i] = $t2; } # A × U for $i (0..N) { $t1 = $$a[$i][$p] * $c + $$a[$i][$q] * $s; $t2 = -$$a[$i][$p] * $s + $$a[$i][$q] * $c; $$a[$i][$p] = $t1; $$a[$i][$q] = $t2; } # 行列の対角要素を表示 (固有値) printf("%3d\t", $k); disp_eigenvalue(\@$a); # 収束判定 last if ($max_val < 0.00000000001); } } # 対角要素を表示 sub disp_eigenvalue { my ($a) = @_; for $i (0..N) { printf("%14.10f\t", $$a[$i][$i]); } print "\n"; } # 固有ベクトルを表示 sub disp_eigenvector { my ($matrix) = @_; for $row (@$matrix) { normarize(\@$row); disp_vector(\@$row); } } # 1次元配列を表示 sub disp_vector { my ($x) = @_; for $i (0..N) { printf("%14.10f\t", $$x[$i]); } print "\n"; } # 正規化 (ベクトルの長さを1にする) sub normarize { my ($x) = @_; my $s = 0.0; for $i (0..N) { $s += $$x[$i] * $$x[$i]; } $s = sqrt($s); for $i (0..N) { $$x[$i] /= $s; } }
Z:\>perl Z:\1105.pl 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
PHP
<?php define("N", 3); # ヤコビ法で固有値を求める main(); function main() { $a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]]; $v = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]]; # ヤコビ法 jacobi($a, $v); print "\neigenvalue\n"; disp_eigenvalue($a); print "\neigenvector\n"; disp_eigenvector($v); } # ヤコビ法 function jacobi(&$a, &$v) { foreach (range(1, 100) as $k) { # 最大値を探す $max_val = 0.0; foreach (range(0, (N - 1)) as $i) { foreach (range(($i + 1), N) as $j) { if ($max_val < abs($a[$i][$j])) { $max_val = abs($a[$i][$j]); $p = $i; $q = $j; } } } # θ を求める if (abs($a[$p][$p] - $a[$q][$q]) < 0.00000000001) { # a_{pp} = a_{qq} のとき、回転角tをπ/4にする $t = M_PI / 4.0; if ($a[$p][$p] < 0.0) $t = -$t; } else { # a_{pp} ≠ a_{qq} のとき $t = atan(2.0 * $a[$p][$q] / ($a[$p][$p] - $a[$q][$q])) / 2.0; } # θ を使って 行列 U を作成し、A = U^t × A × U $c = cos($t); $s = sin($t); # U^t × A foreach (range(0, N) as $i) { $t1 = $a[$p][$i] * $c + $a[$q][$i] * $s; $t2 = -$a[$p][$i] * $s + $a[$q][$i] * $c; $a[$p][$i] = $t1; $a[$q][$i] = $t2; # 固有ベクトル $t1 = $v[$p][$i] * $c + $v[$q][$i] * $s; $t2 = -$v[$p][$i] * $s + $v[$q][$i] * $c; $v[$p][$i] = $t1; $v[$q][$i] = $t2; } # A × U foreach (range(0, N) as $i) { $t1 = $a[$i][$p] * $c + $a[$i][$q] * $s; $t2 = -$a[$i][$p] * $s + $a[$i][$q] * $c; $a[$i][$p] = $t1; $a[$i][$q] = $t2; } # 行列の対角要素を表示 (固有値) printf("%3d\t", $k); disp_eigenvalue($a); # 収束判定 if ($max_val < 0.00000000001) break; } } # 対角要素を表示 function disp_eigenvalue($a) { foreach (range(0, N) as $i) printf("%14.10f\t", $a[$i][$i]); print "\n"; } # 1次元配列を表示 function disp_vector($row) { foreach ($row as $col) printf("%14.10f\t", $col); print "\n"; } # 正規化 (ベクトルの長さを1にする) function normarize(&$x) { $s = 0.0; foreach (range(0, N) as $i) { $s += $x[$i] * $x[$i]; } $s = sqrt($s); foreach (range(0, N) as $i) { $x[$i] /= $s; } } # 固有ベクトルを表示 function disp_eigenvector($matrix) { foreach ($matrix as $row) { normarize($row); disp_vector($row); } } ?>
Z:\>php Z:\1105.php 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
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 "" # 1次元配列を表示 def disp_vector(row): for col in row: print "%14.10f\t" % col, print "" # 正規化 (ベクトルの長さを1にする) def normarize(x): s = 0.0 for i in range(0, N, 1): s += x[i] * x[i] s = math.sqrt(s) for i in range(0, N, 1): x[i] /= s # 固有ベクトルを表示 def disp_eigenvector(matrix): for i in range(0, N, 1): row = [0.0 ,0.0 ,0.0 ,0.0] for j in range(0, N, 1): row[j] = matrix[i][j] normarize(row) disp_vector(row) # ヤコビ法 def jacobi(a, v): for k in range(1, 200, 1): # 最大値を探す p = 0 q = 0 max_val = 0.0 for i in range(0, N, 1): for j in range(i + 1, N, 1): if (max_val < abs(a[i][j])): max_val = abs(a[i][j]) p = i q = j # θ を求める t = 0.0 if (abs(a[p][p] - a[q][q]) < 0.00000000001): # a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = math.pi / 4.0 if (a[p][p] < 0): t = -t else: # a_{pp} ≠ a_{qq} のとき t = math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 # θ を使って 行列 U を作成し、A = U^t × A × U c = math.cos(t) s = math.sin(t) # U^t × A t1 = 0.0 t2 = 0.0 for i in range(0, N, 1): t1 = a[p][i] * c + a[q][i] * s t2 = -a[p][i] * s + a[q][i] * c a[p][i] = t1 a[q][i] = t2 # 固有ベクトル t1 = v[p][i] * c + v[q][i] * s t2 = -v[p][i] * s + v[q][i] * c v[p][i] = t1 v[q][i] = t2 # A × U for i in range(0, N, 1): t1 = a[i][p] * c + a[i][q] * s t2 = -a[i][p] * s + a[i][q] * c a[i][p] = t1 a[i][q] = t2 # 対角要素を表示 print "%3d\t" % k, disp_eigenvalue(a) # 収束判定 if (max_val < 0.00000000001): break # ヤコビ法で固有値を求める 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]] v = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]] # ヤコビ法 jacobi(a, v) print "" print "eigenvalue" disp_eigenvalue(a) print "" print "eigenvector" disp_eigenvector(v)
Z:\>python Z:\1105.py 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
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 # 1次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # 正規化 (ベクトルの長さを1にする) def normarize(x) s = 0.0 (0..N).each do |i| s += x[i] * x[i] end s = Math.sqrt(s) (0..N).each do |i| x[i] /= s end end # 固有ベクトルを表示 def disp_eigenvector(matrix) (0..N).each do |i| row = [0.0 ,0.0 ,0.0 ,0.0] (0..N).each do |j| row[j] = matrix[i][j] end normarize(row) disp_vector(row) end end # ヤコビ法 def jacobi(a, v) (1..200).each do |k| # 最大値を探す p = 0 q = 0 max_val = 0.0 (0..N).each do |i| ((i + 1)..N).each do |j| if (max_val < (a[i][j]).abs) max_val = (a[i][j]).abs p = i q = j end end end # θ を求める t = 0.0 if ((a[p][p] - a[q][q]).abs < 0.00000000001) # a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math::PI / 4.0 t = -t if (a[p][p] < 0) else # a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 end # θ を使って 行列 U を作成し、A = U^t × A × U c = Math.cos(t) s = Math.sin(t) # U^t × A t1 = 0.0 t2 = 0.0 (0..N).each do |i| t1 = a[p][i] * c + a[q][i] * s t2 = -a[p][i] * s + a[q][i] * c a[p][i] = t1 a[q][i] = t2 # 固有ベクトル t1 = v[p][i] * c + v[q][i] * s t2 = -v[p][i] * s + v[q][i] * c v[p][i] = t1 v[q][i] = t2 end # A × U (0..N).each do |i| t1 = a[i][p] * c + a[i][q] * s t2 = -a[i][p] * s + a[i][q] * c a[i][p] = t1 a[i][q] = t2 end # 対角要素を表示 printf("%3d\t", k) disp_eigenvalue(a) # 収束判定 break if (max_val < 0.00000000001) end end # ヤコビ法で固有値を求める a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] v = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]] # ヤコビ法 jacobi(a, v) puts "" puts "eigenvalue" disp_eigenvalue(a) puts "" puts "eigenvector" disp_eigenvector(v)
Z:\>ruby Z:\1105.rb 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
Groovy
N = 3 // ヤコビ法で固有値を求める main() def main() { def a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] as double[][] def v = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]] as double[][] // ヤコビ法 jacobi(a, v) println() println("eigenvalue") disp_eigenvalue(a) println() println("eigenvector") disp_eigenvector(v) } // ヤコビ法 def jacobi(a, v) { for (k in 1..200) { // 最大値を探す p = 0 q = 0 max_val = 0.0 for (i in 0..(N - 1)) { for (j in (i + 1)..N) { if (max_val < Math.abs(a[i][j])) { max_val = Math.abs(a[i][j]) p = i q = j } } } // θ を求める t = 0.0 if (Math.abs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 if (a[p][p] < 0) t = -t } else { // a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 } // θ を使って 行列 U を作成し、A = U^t × A × U c = Math.cos(t) s = Math.sin(t) // U^t × A t1 = 0.0 t2 = 0.0 for (i in 0..N) { t1 = a[p][i] * c + a[q][i] * s t2 = -a[p][i] * s + a[q][i] * c a[p][i] = t1 a[q][i] = t2 // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s t2 = -v[p][i] * s + v[q][i] * c v[p][i] = t1 v[q][i] = t2 } // A × U for (i in 0..N) { t1 = a[i][p] * c + a[i][q] * s t2 = -a[i][p] * s + a[i][q] * c a[i][p] = t1 a[i][q] = t2 } // 対角要素を表示 printf("%3d\t" , k) disp_eigenvalue(a) // 収束判定 if (max_val < 0.00000000001) break } } // 対角要素を表示 def disp_eigenvalue(a) { for (i in 0..N) printf("%14.10f\t" , a[i][i]) println() } // 固有ベクトルを表示 def disp_eigenvector(matrix) { for (i in 0..N) { row = [0.0 ,0.0 ,0.0 ,0.0] for (j in 0..N) row[j] = matrix[i][j] normarize(row) disp_vector(row) } } // 1次元配列を表示 def disp_vector(row) { for (col in row) { printf("%14.10f\t" , col) } println() } // 正規化 (ベクトルの長さを1にする) def normarize(x) { s = 0.0 for (i in 0..N) s += x[i] * x[i] s = Math.sqrt(s) for (i in 0..N) x[i] /= s }
Z:\>gvy Gvy1105.gvy 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
Pascal
program Pas1105(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; // 正規化 (ベクトルの長さを1にする) procedure normarize(var x:array of Double); var s:Double = 0.0; i:Integer; begin for i := Low(x) to High(x) do s := s + x[i] * x[i]; s := Sqrt(s); for i := Low(x) to High(x) do x[i] := x[i] / s; end; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 固有ベクトルを表示 procedure disp_eigenvector(matrix:TwoDimArray); var i, j:Integer; row:array [0..N] of Double; begin for i := 0 to N do begin for j := 0 to N do row[j] := matrix[i][j]; normarize(row); disp_vector(row); end; end; // ヤコビ法 procedure jacobi(var a:TwoDimArray; var v:TwoDimArray); var i, j, k:Integer; p, q:Integer; max_val:Double; c, s:Double; t, t1, t2:Double; begin for k := 1 to 100 do begin // 最大値を探す max_val := 0.0; for i := 0 to N do begin for j := i + 1 to N do begin if (max_val < Abs(a[i][j])) then begin max_val := Abs(a[i][j]); p := i; q := j; end; end; end; // θ を求める t := 0.0; if Abs(a[p][p] - a[q][q]) < 0.00000000001 then begin // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t := PI / 4.0; if (a[p][p] < 0) then t := -t; end else begin // a_{pp} ≠ a_{qq} のとき t := Arctan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0; end; // θ を使って 行列 U を作成し、A = U^t × A × U c := Cos(t); s := Sin(t); // U^t × A t1 := 0.0; t2 := 0.0; for i := 0 to N do begin t1 := a[p][i] * c + a[q][i] * s; t2 := -a[p][i] * s + a[q][i] * c; a[p][i] := t1; a[q][i] := t2; // 固有ベクトル t1 := v[p][i] * c + v[q][i] * s; t2 := -v[p][i] * s + v[q][i] * c; v[p][i] := t1; v[q][i] := t2; end; // A × U for i := 0 to N do begin t1 := a[i][p] * c + a[i][q] * s; t2 := -a[i][p] * s + a[i][q] * c; a[i][p] := t1; a[i][q] := t2; end; // 対角要素を表示 write(format('%3d'#9, [k])); disp_eigenvalue(a); // 収束判定 if max_val < 0.00000000001 then break; end; end; // ヤコビ法で固有値を求める var a :TwoDimArray = ((5.0, 4.0, 1.0, 1.0) ,(4.0, 5.0, 1.0, 1.0) ,(1.0, 1.0, 4.0, 2.0) ,(1.0, 1.0, 2.0, 4.0)); v :TwoDimArray = ((1.0 ,0.0 ,0.0 ,0.0), (0.0 ,1.0 ,0.0 ,0.0), (0.0 ,0.0 ,1.0 ,0.0), (0.0 ,0.0 ,0.0 ,1.0)); begin // ヤコビ法 jacobi(a, v); writeln(); writeln('eigenvalue'); disp_eigenvalue(a); writeln(); writeln('eigenvector'); disp_eigenvector(v); end.
Z:\>fpc -v0 -l- Pas1105.pp Z:\>Pas1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
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 Ada1105 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)); v:Long_Float_TwoDimArray := ((1.0 ,0.0 ,0.0 ,0.0) ,(0.0 ,1.0 ,0.0 ,0.0) ,(0.0 ,0.0 ,1.0 ,0.0) ,(0.0 ,0.0 ,0.0 ,1.0)); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 正規化 (ベクトルの長さを1にする) procedure normarize(x:in out Long_Float_Array) is s:Long_Float; begin s := 0.0; for i in x'Range loop s := s + x(i) * x(i); end loop; s := Sqrt(s); for i in x'Range loop x(i) := x(i) / s; end loop; end normarize; -- 対角要素を表示 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; New_Line; end disp_eigenvalue; -- 固有ベクトルを表示 procedure disp_eigenvector(matrix:Long_Float_TwoDimArray) is row:Long_Float_Array; begin for i in 0..N loop for j in 0..N loop row(j) := matrix(i, j); end loop; normarize(row); disp_vector(row); end loop; end disp_eigenvector; -- ヤコビ法 procedure jacobi(a:in out Long_Float_TwoDimArray; v:in out Long_Float_TwoDimArray) is p, q:Integer; max_val:Long_Float; c, s, t, t1, t2:Long_Float; begin for k in 1..100 loop -- 最大値を探す max_val := 0.0; for i in 0..N loop for j in (i + 1)..N loop if (max_val < Abs(a(i, j))) then max_val := Abs(a(i, j)); p := i; q := j; end if; end loop; end loop; -- θ を求める t := 0.0; if Abs(a(p, p) - a(q, q)) < 0.00000000001 then -- a_{pp} = a_{qq} のとき、回転角tをπ/4にする t := Ada.Numerics.Pi / 4.0; if (a(p, p) < 0.0) then t := -t; end if; else -- a_{pp} ≠ a_{qq} のとき t := Arctan(2.0 * a(p, q) / (a(p, p) - a(q, q))) / 2.0; end if; -- θ を使って 行列 U を作成し、A = U^t × A × U c := Cos(t); s := Sin(t); -- U^t × A t1 := 0.0; t2 := 0.0; for i in 0..N loop t1 := a(p, i) * c + a(q, i) * s; t2 := -a(p, i) * s + a(q, i) * c; a(p, i) := t1; a(q, i) := t2; -- 固有ベクトル t1 := v(p, i) * c + v(q, i) * s; t2 := -v(p, i) * s + v(q, i) * c; v(p, i) := t1; v(q, i) := t2; end loop; -- A × U for i in 0..N loop t1 := a(i, p) * c + a(i, q) * s; t2 := -a(i, p) * s + a(i, q) * c; a(i, p) := t1; a(i, q) := t2; end loop; --対角要素を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_eigenvalue(a); -- 収束判定 if max_val < 0.00000000001 then exit; end if; end loop; end jacobi; -- ヤコビ法で固有値を求める begin -- ヤコビ法 jacobi(a, v); Put_Line(""); Put_Line("eigenvalue"); disp_eigenvalue(a); New_Line; Put_Line("eigenvector"); disp_eigenvector(v); end Ada1105;
xxxxxx@yyyyyy /Z $ gnatmake Ada1105.adb xxxxxx@yyyyyy /Z $ Ada1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
VB.NET
Option Explicit Module VB1105 Private Const N As Integer = 3 'ヤコビ法で固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} Dim v(,) As Double = {{1.0 ,0.0 ,0.0 ,0.0}, {0.0 ,1.0 ,0.0 ,0.0}, {0.0 ,0.0 ,1.0 ,0.0}, {0.0 ,0.0 ,0.0 ,1.0}} 'ヤコビ法 jacobi(a, v) Console.WriteLine() Console.WriteLine("eigenvalue") disp_eigenvalue(a) Console.WriteLine() Console.WriteLine("eigenvector") disp_eigenvector(v) End Sub 'ヤコビ法 Private Sub jacobi(ByRef a(,) As Double, ByRef v(,) As Double) For k As Integer = 1 To 100 '最大値を探す Dim p, q As Integer Dim max_val As Double = 0.0 For i As Integer = 0 To N For j As Integer = i + 1 To N If (max_val < Math.Abs(a(i, j))) Then max_val = Math.Abs(a(i, j)) p = i q = j End If Next Next 'θ を求める Dim t As Double = 0.0 If Math.Abs(a(p, p) - a(q, q)) < 0.00000000001 Then 'a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 If (a(p, p) < 0) Then t = -t End If Else 'a_{pp} ≠ a_{qq} のとき t = Math.Atan(2.0 * a(p, q) / (a(p, p) - a(q, q))) / 2.0 End If 'θ を使って 行列 U を作成し、A = U^t × A × U Dim c As Double = Math.Cos(t) Dim s As Double = Math.Sin(t) 'U^t × A Dim t1 As Double = 0.0 Dim t2 As Double = 0.0 For i As Integer = 0 To N t1 = a(p, i) * c + a(q, i) * s t2 = -a(p, i) * s + a(q, i) * c a(p, i) = t1 a(q, i) = t2 '固有ベクトル t1 = v(p, i) * c + v(q, i) * s t2 = -v(p, i) * s + v(q, i) * c v(p, i) = t1 v(q, i) = t2 Next 'A × U For i As Integer = 0 To N t1 = a(i, p) * c + a(i, q) * s t2 = -a(i, p) * s + a(i, q) * c a(i, p) = t1 a(i, q) = t2 Next '対角要素を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_eigenvalue(a) '収束判定 If max_val < 0.00000000001 Then Exit For 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 '固有ベクトルを表示 Private Sub disp_eigenvector(ByVal matrix(,) As Double) For i As Integer = 0 To N Dim row() As Double = New Double(N){} For j As Integer = 0 To N row(j) = matrix(i, j) Next normarize(row) disp_vector(row) Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x() As Double) Dim s As Double = 0.0 For i As Integer = 0 To N s += x(i) * x(i) Next s = Math.Sqrt(s) For i As Integer = 0 To N x(i) /= s Next End Sub End Module
Z:\>vbc -nologo VB1105.vb Z:\>VB1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
C#
using System; public class CS1105 { private const int N = 4; // ヤコビ法で固有値を求める public static void Main() { double[,] a = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double[,] v = {{1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 1.0}}; // ヤコビ法 jacobi(a, v); Console.WriteLine(); Console.WriteLine("eigenvalue"); disp_eigenvalue(a); Console.WriteLine(); Console.WriteLine("eigenvector"); disp_eigenvector(v); } // ヤコビ法 private static void jacobi(double[,] a, double[,] v) { for (int k = 1; k <= 100; k++) { // 最大値を探す int p = 0; int q = 0; double max_val = 0.0; for (int i = 0; i < N; i++) { for (int j = i + 1; j < N; j++) { if (max_val < Math.Abs(a[i, j])) { max_val = Math.Abs(a[i, j]); p = i; q = j; } } } // θ を求める double t = 0.0; if (Math.Abs(a[p, p] - a[q, q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0; if (a[p, p] < 0) t = -t; } else { // a_{pp} ≠ a_{qq} のとき t = Math.Atan(2.0 * a[p, q] / (a[p, p] - a[q, q])) / 2.0; } // θ を使って 行列 U を作成し、A = U^t × A × U double c = Math.Cos(t); double s = Math.Sin(t); // U^t × A double t1 = 0.0; double t2 = 0.0; for (int i = 0; i < N; i++) { t1 = a[p, i] * c + a[q, i] * s; t2 = -a[p, i] * s + a[q, i] * c; a[p, i] = t1; a[q, i] = t2; // 固有ベクトル t1 = v[p, i] * c + v[q, i] * s; t2 = -v[p, i] * s + v[q, i] * c; v[p, i] = t1; v[q, i] = t2; } // A × U for (int i = 0; i < N; i++) { t1 = a[i, p] * c + a[i, q] * s; t2 = -a[i, p] * s + a[i, q] * c; a[i, p] = t1; a[i, q] = t2; } // 対角要素を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_eigenvalue(a); // 収束判定 if (max_val < 0.00000000001) break; } } // 対角要素を表示 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(); } // 固有ベクトルを表示 private static void disp_eigenvector(double[,] matrix) { for (int i = 0; i < N; i++) { double[] row = new double[N]; for (int j = 0; j < N; j++) row[j] = matrix[i, j] ; normarize(row); disp_vector(row); } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 正規化 (ベクトルの長さを1にする) private static void normarize(double[] x) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = Math.Sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; } }
Z:\>csc -nologo CS1105.cs Z:\>CS1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
Java
import java.lang.*; public class Java1105 { private static final int N = 4; // ヤコビ法で固有値を求める public static void main(String []args) { double[][] a = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double[][] v = {{1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 1.0}}; // ヤコビ法 jacobi(a, v); System.out.println(); System.out.println("eigenvalue"); disp_eigenvalue(a); System.out.println(); System.out.println("eigenvector"); disp_eigenvector(v); } // ヤコビ法 private static void jacobi(double[][] a, double[][] v) { for (int k = 1; k <= 100; k++) { // 最大値を探す int p = 0; int q = 0; double max_val = 0.0; for (int i = 0; i < N; i++) { for (int j = i + 1; j < N; j++) { if (max_val < Math.abs(a[i][j])) { max_val = Math.abs(a[i][j]); p = i; q = j; } } } // θ を求める double t = 0.0; if (Math.abs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0; if (a[p][p] < 0) t = -t; } else { // a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0; } // θ を使って 行列 U を作成し、A = U^t × A × U double c = Math.cos(t); double s = Math.sin(t); // U^t × A double t1 = 0.0; double t2 = 0.0; for (int i = 0; i < N; i++) { t1 = a[p][i] * c + a[q][i] * s; t2 = -a[p][i] * s + a[q][i] * c; a[p][i] = t1; a[q][i] = t2; // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s; t2 = -v[p][i] * s + v[q][i] * c; v[p][i] = t1; v[q][i] = t2; } // A × U for (int i = 0; i < N; i++) { t1 = a[i][p] * c + a[i][q] * s; t2 = -a[i][p] * s + a[i][q] * c; a[i][p] = t1; a[i][q] = t2; } // 対角要素を表示 System.out.print(String.format("%3d\t", k)); disp_eigenvalue(a); // 収束判定 if (max_val < 0.00000000001) break; } } // 対角要素を表示 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(); } // 固有ベクトルを表示 private static void disp_eigenvector(double[][] matrix) { for (int i = 0; i < N; i++) { double[] row = new double[N]; for (int j = 0; j < N; j++) row[j] = matrix[i][j] ; normarize(row); disp_vector(row); } } // 1次元配列を表示 private static void disp_vector(double[] row) { for (double col: row) System.out.print(String.format("%14.10f\t", col)); System.out.println(); } // 正規化 (ベクトルの長さを1にする) private static void normarize(double[] x) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = Math.sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; } }
Z:\>javac Java1105.java Z:\>java Java1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
C++
#include <iostream> #include <iomanip> #include <math.h> using namespace std; const int N = 4; // ヤコビ法 void jacobi(double a[N][N], double v[N][N]); // 対角要素を表示 void disp_eigenvalue(double a[N][N]); // 固有ベクトルを表示 void disp_eigenvector(double matrix[N][N]); // 1次元配列を表示 void disp_vector(double row[N]); // 正規化 (ベクトルの長さを1にする) void normarize(double x[N]); // ヤコビ法で固有値を求める int main() { double a[N][N] = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double v[N][N] = {{1.0 ,0.0 ,0.0 ,0.0}, {0.0 ,1.0 ,0.0 ,0.0}, {0.0 ,0.0 ,1.0 ,0.0}, {0.0 ,0.0 ,0.0 ,1.0}}; // ヤコビ法 jacobi(a, v); cout << endl << "eigenvalue" << endl; disp_eigenvalue(a); cout << endl << "eigenvector" << endl; disp_eigenvector(v); return 0; } // ヤコビ法 void jacobi(double a[N][N], double v[N][N]) { for (int k = 1; k <= 100; k++) { // 最大値を探す int p = 0; int q = 0; double max_val = 0.0; for (int i = 0; i < N; i++) { for (int j = i + 1; j < N; j++) { if (max_val < fabs(a[i][j])) { max_val = fabs(a[i][j]); p = i; q = j; } } } // θ を求める double t; if (fabs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = M_PI / 4.0; if (a[p][p] < 0) t = -t; } else { // a_{pp} ≠ a_{qq} のとき t = atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0; } // θ を使って 行列 U を作成し、A = U^t × A × U double c = cos(t); double s = sin(t); // U^t × A double t1; double t2; for (int i = 0; i < N; i++) { t1 = a[p][i] * c + a[q][i] * s; t2 = -a[p][i] * s + a[q][i] * c; a[p][i] = t1; a[q][i] = t2; // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s; t2 = -v[p][i] * s + v[q][i] * c; v[p][i] = t1; v[q][i] = t2; } // A × U for (int i = 0; i < N; i++) { t1 = a[i][p] * c + a[i][q] * s; t2 = -a[i][p] * s + a[i][q] * c; a[i][p] = t1; a[i][q] = t2; } // 対角要素を表示 cout << setw(3) << k << "\t"; disp_eigenvalue(a); // 収束判定 if (max_val < 0.00000000001) break; } } // 対角要素を表示 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; } // 固有ベクトルを表示 void disp_eigenvector(double matrix[N][N]) { for (int i = 0; i < N; i++) { double row[N]; for (int j = 0; j < N; j++) row[j] = matrix[i][j] ; normarize(row); disp_vector(row); } } // 1次元配列を表示 void disp_vector(double row[N]) { for (int i = 0; i < N; i++) cout << setw(14) << fixed << setprecision(10) << row[i] << "\t"; cout << endl; } // 正規化 (ベクトルの長さを1にする) void normarize(double x[N]) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; }
Z:\>bcc32 -q CP1105.cpp cp1105.cpp: Z:\>CP1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
Objective-C
#import <Foundation/Foundation.h> #import <math.h> #define N 4 // ヤコビ法 void jacobi(double a[N][N], double v[N][N]); // 対角要素を表示 void disp_eigenvalue(double a[N][N]); // 固有ベクトルを表示 void disp_eigenvector(double matrix[N][N]); // 1次元配列を表示 void disp_vector(double row[N]); // 正規化 (ベクトルの長さを1にする) void normarize(double x[N]); // ヤコビ法で固有値を求める int main() { double a[N][N] = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}}; double v[N][N] = {{1.0 ,0.0 ,0.0 ,0.0}, {0.0 ,1.0 ,0.0 ,0.0}, {0.0 ,0.0 ,1.0 ,0.0}, {0.0 ,0.0 ,0.0 ,1.0}}; // ヤコビ法 jacobi(a, v); printf("\neigenvalue\n"); disp_eigenvalue(a); printf("\neigenvector\n"); disp_eigenvector(v); return 0; } // ヤコビ法 void jacobi(double a[N][N], double v[N][N]) { int i, j, k; for (k = 1; k <= 100; k++) { // 最大値を探す int p = 0; int q = 0; double max_val = 0.0; for (i = 0; i < N; i++) { for (j = i + 1; j < N; j++) { if (max_val < fabs(a[i][j])) { max_val = fabs(a[i][j]); p = i; q = j; } } } // θ を求める double t = 0.0; if (fabs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = M_PI / 4.0; if (a[p][p] < 0) t = -t; } else { // a_{pp} ≠ a_{qq} のとき t = atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0; } // θ を使って 行列 U を作成し、A = U^t × A × U double c = cos(t); double s = sin(t); // U^t × A double t1 = 0.0; double t2 = 0.0; for (i = 0; i < N; i++) { t1 = a[p][i] * c + a[q][i] * s; t2 = -a[p][i] * s + a[q][i] * c; a[p][i] = t1; a[q][i] = t2; // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s; t2 = -v[p][i] * s + v[q][i] * c; v[p][i] = t1; v[q][i] = t2; } // A × U for (i = 0; i < N; i++) { t1 = a[i][p] * c + a[i][q] * s; t2 = -a[i][p] * s + a[i][q] * c; a[i][p] = t1; a[i][q] = t2; } // 対角要素を表示 printf("%3d\t", k); disp_eigenvalue(a); // 収束判定 if (max_val < 0.00000000001) break; } } // 対角要素を表示 void disp_eigenvalue(double a[N][N]) { int i; for (i = 0; i < N; i++) printf("%14.10f\t", a[i][i]); printf("\n"); } // 固有ベクトルを表示 void disp_eigenvector(double matrix[N][N]) { int i, j; for (i = 0; i < N; i++) { double row[N]; for (j = 0; j < N; j++) row[j] = matrix[i][j] ; normarize(row); disp_vector(row); } } // 1次元配列を表示 void disp_vector(double row[N]) { int i; for (i = 0; i < N; i++) printf("%14.10f\t", row[i]); printf("\n"); } // 正規化 (ベクトルの長さを1にする) void normarize(double x[N]) { int i; double s = 0.0; for (i = 0; i < N; i++) s += x[i] * x[i]; s = sqrt(s); for (i = 0; i < N; i++) x[i] /= s; }
xxxxxx@yyyyyy /Z $ gcc -o OC1105 OC1105.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
D
import std.stdio; import std.math; const int N = 4; // ヤコビ法で固有値を求める void main(string[] args) { double a[N][N] = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]]; double v[N][N] = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]]; // ヤコビ法 jacobi(a, v); writefln("\neigenvalue"); disp_eigenvalue(a); writefln("\neigenvector"); disp_eigenvector(v); } // ヤコビ法 void jacobi(ref double[N][N] a, ref double[N][N] v) { foreach (k; 1..100) { // 最大値を探す int p = 0; int q = 0; double max_val = 0.0; foreach (i; 0..N) { foreach (j; i + 1..N) { if (max_val < fabs(a[i][j])) { max_val = fabs(a[i][j]); p = i; q = j; } } } // θ を求める double t = 0.0; if (fabs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = PI / 4.0; if (a[p][p] < 0) t = -t; } else { // a_{pp} ≠ a_{qq} のとき t = atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0; } // θ を使って 行列 U を作成し、A = U^t × A × U double c = cos(t); double s = sin(t); // U^t × A double t1 = 0.0; double t2 = 0.0; foreach (i; 0..N) { t1 = a[p][i] * c + a[q][i] * s; t2 = -a[p][i] * s + a[q][i] * c; a[p][i] = t1; a[q][i] = t2; // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s; t2 = -v[p][i] * s + v[q][i] * c; v[p][i] = t1; v[q][i] = t2; } // A × U foreach (i; 0..N) { t1 = a[i][p] * c + a[i][q] * s; t2 = -a[i][p] * s + a[i][q] * c; a[i][p] = t1; a[i][q] = t2; } // 対角要素を表示 writef("%3d\t", k); disp_eigenvalue(a); // 収束判定 if (max_val < 0.00000000001) break; } } // 対角要素を表示 void disp_eigenvalue(double[N][N] a) { foreach (i; 0..N) writef("%14.10f\t", a[i][i]); writefln(""); } // 固有ベクトルを表示 void disp_eigenvector(double[N][N] matrix) { foreach (i; 0..N) { double row[N]; foreach (j; 0..N) row[j] = matrix[i][j] ; normarize(row); disp_vector(row); } } // 1次元配列を表示 void disp_vector(double[N] row) { foreach (i; 0..N) writef("%14.10f\t", row[i]); writefln(""); } // 正規化 (ベクトルの長さを1にする) void normarize(ref double[N] x) { double s = 0.0; foreach (i; 0..N) s += x[i] * x[i]; s = sqrt(s); foreach (i; 0..N) x[i] /= s; }
Z:\>dmd D1105.d Z:\>D1105 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
Go
package main import "fmt" import "math" const N = 4 // ヤコビ法で固有値を求める func main() { var a [N][N]float64 = [N][N]float64{{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} var v [N][N]float64 = [N][N]float64{{1.0 ,0.0 ,0.0 ,0.0}, {0.0 ,1.0 ,0.0 ,0.0}, {0.0 ,0.0 ,1.0 ,0.0}, {0.0 ,0.0 ,0.0 ,1.0}} // ヤコビ法 jacobi(&a, &v) fmt.Println("\neigenvalue") disp_eigenvalue(a) fmt.Println("\neigenvector") disp_eigenvector(v) } // ヤコビ法 func jacobi(a *[N][N]float64, v *[N][N]float64) { for k := 1; k < 100; k++ { // 最大値を探す var p int = 0 var q int = 0 var max_val float64 = 0.0 for i := 0; i < N; i++ { for j := i + 1; j < N; j++ { if (max_val < math.Abs(a[i][j])) { max_val = math.Abs(a[i][j]) p = i q = j } } } // θ を求める var t float64 = 0.0 if (math.Abs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角θをπ/4にする t = math.Pi / 4.0 if (a[p][p] < 0) { t = -t } } else { // a_{pp} ≠ a_{qq} のとき t = math.Atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 } // θ を使って 行列 P を作成し、A = P^t × A × P var c float64 = math.Cos(t) var s float64 = math.Sin(t) // P^t × A var t1 float64 = 0.0 var t2 float64 = 0.0 for i := 0; i < N; i++ { t1 = a[p][i] * c + a[q][i] * s t2 = -a[p][i] * s + a[q][i] * c a[p][i] = t1 a[q][i] = t2 // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s t2 = -v[p][i] * s + v[q][i] * c v[p][i] = t1 v[q][i] = t2 } // A × P for i := 0; i < N; i++ { t1 = a[i][p] * c + a[i][q] * s t2 = -a[i][p] * s + a[i][q] * c a[i][p] = t1 a[i][q] = t2 } // 対角要素を表示 fmt.Printf("%3d\t", k) disp_eigenvalue(*a) // 収束判定 if (max_val < 0.00000000001) { break } } } // 対角要素を表示 func disp_eigenvalue(a[N][N]float64) { for i := 0; i < N; i++ { fmt.Printf("%14.10f\t", a[i][i]) } fmt.Println("") } // 固有ベクトルを表示 func disp_eigenvector(matrix[N][N]float64) { for i := 0; i < N; i++ { var row = []float64 {0.0 ,0.0 ,0.0 ,0.0} for j := 0; j < N; j++ { row[j] = matrix[i][j] } normarize(row) disp_vector(row) } } // 1次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // 正規化 (ベクトルの長さを1にする) func normarize(x[]float64) { var s float64 = 0.0 for i := 0; i < N; i++ { s += x[i] * x[i] } s = math.Sqrt(s) for i := 0; i < N; i++ { x[i] /= s } }
Z:\>go run GO1105.go 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
Scala
object Scala1105 { val N = 3 // ヤコビ法で固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) var v:Array[Array[Double]] = Array(Array(1.0, 0.0, 0.0, 0.0) ,Array(0.0, 1.0, 0.0, 0.0) ,Array(0.0, 0.0, 1.0, 0.0) ,Array(0.0, 0.0, 0.0, 1.0)) // ヤコビ法 val (a1, v1) = jacobi(a, v) println() println("eigenvalue") disp_eigenvalue(a1) println() println("eigenvector") disp_eigenvector(v1) } // ヤコビ法 def jacobi(a: => Array[Array[Double]], v: => Array[Array[Double]]): (Array[Array[Double]], Array[Array[Double]]) = { var max_val = 0.0 var k = 1 do { // 最大値を探す var p = 0 var q = 0 max_val = 0.0 for (i <- 0 to N) { for (j <- i + 1 to N) { if (max_val < Math.abs(a(i)(j))) { max_val = Math.abs(a(i)(j)) p = i q = j } } } // θ を求める var t = 0.0 if (Math.abs(a(p)(p) - a(q)(q)) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 if (a(p)(p) < 0) t = -t } else { // a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a(p)(q) / (a(p)(p) - a(q)(q))) / 2.0 } // θ を使って 行列 U を作成し、A = U^t × A × U val c = Math.cos(t) val s = Math.sin(t) // U^t × A var t1 = 0.0 var t2 = 0.0 for (i <- 0 to N) { t1 = a(p)(i) * c + a(q)(i) * s t2 = -a(p)(i) * s + a(q)(i) * c a(p)(i) = t1 a(q)(i) = t2 // 固有ベクトル t1 = v(p)(i) * c + v(q)(i) * s t2 = -v(p)(i) * s + v(q)(i) * c v(p)(i) = t1 v(q)(i) = t2 } // A × U for (i <- 0 to N) { t1 = a(i)(p) * c + a(i)(q) * s t2 = -a(i)(p) * s + a(i)(q) * c a(i)(p) = t1 a(i)(q) = t2 } // 対角要素を表示 print("%3d\t".format(k)) disp_eigenvalue(a) k += 1 } while (k <= 100 && max_val >= 0.00000000001) return (a, v) } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.map( col => print("%14.10f\t".format(col))) println() } // 正規化 (ベクトルの長さを1にする) def normarize(x:Array[Double]):Array[Double] = { x.map(n => n / Math.sqrt(x.map(n => n * n).sum)) } // 対角要素を表示 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() } // 固有ベクトルを表示 def disp_eigenvector(matrix:Array[Array[Double]]) { for (i <- 0 to N) { disp_vector(normarize((for (j <- 0 to N) yield (matrix(i)(j))).toArray)) } } }
Z:\>scala Scala1105.scala 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 -0.0000000000 -0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812
F#
module Fs1105 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 正規化 (ベクトルの長さを1にする) let normarize (x:float[]):float[] = let s = x |> Array.map (fun n -> n * n) |> Array.sum x |> Array.map(fun n -> n / Math.Sqrt(s)) // 対角要素を表示 let disp_eigenvalue (matrix:float[][]) = [| for i in 0..N -> matrix.[i].[i] |] |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 固有ベクトルを表示 let disp_eigenvector (matrix:float[][]) = for i in [0..N] do disp_vector ( normarize [| for j in 0..N -> matrix.[i].[j] |] ) // ヤコビ法 let jacobi (a:float[][]) (v:float[][]) : float[][] * float[][] = let mutable max_val = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 100 && max_val >= 0.00000000001) do if sw then sw <- false // 最大値を探す let mutable p = 0 let mutable q = 0 max_val <- 0.0 for i in [0..N] do for j in [(i + 1)..N] do if max_val < Math.Abs(a.[i].[j]) then max_val <- Math.Abs(a.[i].[j]) p <- i q <- j // θ を求める let mutable t = 0.0 if Math.Abs(a.[p].[p] - a.[q].[q]) < 0.00000000001 then // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t <- Math.PI / 4.0 if a.[p].[p] < 0.0 then t <- -t else // a_{pp} ≠ a_{qq} のとき t <- Math.Atan(2.0 * a.[p].[q] / (a.[p].[p] - a.[q].[q])) / 2.0 // θ を使って 行列 U を作成し、A = U^t × A × U let c = Math.Cos(t) let s = Math.Sin(t) // U^t × A let mutable t1 = 0.0 let mutable t2 = 0.0 for i in [0..N] do t1 <- a.[p].[i] * c + a.[q].[i] * s t2 <- -a.[p].[i] * s + a.[q].[i] * c a.[p].[i] <- t1 a.[q].[i] <- t2 // 固有ベクトル t1 <- v.[p].[i] * c + v.[q].[i] * s t2 <- -v.[p].[i] * s + v.[q].[i] * c v.[p].[i] <- t1 v.[q].[i] <- t2 // A × U for i in [0..N] do t1 <- a.[i].[p] * c + a.[i].[q] * s t2 <- -a.[i].[p] * s + a.[i].[q] * c a.[i].[p] <- t1 a.[i].[q] <- t2 // 対角要素を表示 printf "%3d\t" k disp_eigenvalue a k <- k + 1 (a, v) // ヤコビ法で固有値を求める 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 v:float[][] = [| [|1.0; 0.0; 0.0; 0.0|]; [|0.0; 1.0; 0.0; 0.0|]; [|0.0; 0.0; 1.0; 0.0|]; [|0.0; 0.0; 0.0; 1.0|] |] // ヤコビ法 let (a1, v1) = jacobi a v printfn "" printfn "eigenvalue" disp_eigenvalue a1 printfn "" printfn "eigenvector" disp_eigenvector v1 exit 0
Z:\>fsi --nologo --quiet Fs1105.fs 1 9.0000000000 1.0000000000 4.0000000000 4.0000000000 2 9.0000000000 1.0000000000 6.0000000000 2.0000000000 3 10.0000000000 1.0000000000 5.0000000000 2.0000000000 4 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 -0.7071067812 0.7071067812 0.0000000000 0.0000000000 -0.3162277660 -0.3162277660 0.6324555320 0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812