さまざまな言語で数値計算
Only Do What Only You Can Do
ハウスホルダー法
ハウスホルダー法を使って, 全ての固有値と, 対応する固有ベクトルを求める
対称行列 $ \boldsymbol{A} $ を, 以下のように分割して考える.
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 d: d = Array(0.0 ,0.0 ,0.0 ,0.0) Dim e: e = Array(0.0 ,0.0 ,0.0 ,0.0) 'ハウスホルダー変換 Call tridiagonalize(a, d, e) 'QR分解 Call decomp(a, d, e) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" Call disp_vector(d) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvector" Call disp_matrix(a) End Sub 'ハウスホルダー変換 Private Sub tridiagonalize(ByRef a, ByRef d, ByRef e) Dim v: v = Array(0.0 ,0.0 ,0.0 ,0.0) Dim i, j, k For k = 0 To N - 2 d(k) = a(k)(k) Dim t: t = 0.0 Dim w: w = Array(0.0, 0.0, 0.0, 0.0) For i = k + 1 To N w(i) = a(i)(k) t = t + w(i) * w(i) Next Dim s: s = Sqr(t) If w(k + 1) < 0 Then s = -s End If If Abs(s) < 0.00000000001 Then e(k + 1) = 0.0 Else e(k + 1) = -s w(k + 1) = w(k + 1) + s s = 1 / Sqr(w(k + 1) * s) For i = k + 1 To N w(i) = w(i) * s Next For i = k + 1 To N s = 0.0 For j = k + 1 To N If j <= i Then s = s + a(i)(j) * w(j) Else s = s + a(j)(i) * w(j) End If Next v(i) = s Next s = 0 For i = k + 1 To N s = s + w(i) * v(i) Next s = s / 2 For i = k + 1 To N v(i) = v(i) - s * w(i) Next For i = k + 1 To N For j = k + 1 To i a(i)(j) = a(i)(j) - w(i) * v(j) - w(j) * v(i) Next Next '次の行は固有ベクトルを求めないなら不要 For i = k + 1 To N a(i)(k) = w(i) Next End If Next d(N - 1) = a(N - 1)(N - 1) d(N) = a(N)(N) e(0) = 0.0 e(N) = a(N)(N - 1) '次の行は固有ベクトルを求めないなら不要 For k = N To 0 Step -1 If k < N - 1 Then For i = k + 1 To N w(i) = a(i)(k) Next For i = k + 1 To N s = 0.0 For j = k + 1 To N s = s + a(i)(j) * w(j) Next v(i) = s Next For i = k + 1 To N For j = k + 1 To N a(i)(j) = a(i)(j) - v(i) * w(j) Next Next End If For i = 0 To N a(i)(k) = 0.0 Next a(k)(k) = 1.0 Next End Sub 'QR分解 Private Sub decomp(ByRef a, ByRef d, ByRef e) Dim i, j, k e(0) = 1.0 Dim h: h = N Do While (Abs(e(h)) < 0.00000000001) h = h - 1 Loop Do While h > 0 e(0) = 0.0 Dim l: l = h - 1 Do While (Abs(e(l)) >= 0.00000000001) l = l - 1 Loop For j = 1 To 100 Dim w: w = (d(h - 1) - d(h)) / 2.0 Dim s: s = Sqr(w * w + e(h) * e(h)) If w < 0.0 Then s = - s End If Dim x: x = d(l) - d(h) + e(h) * e(h) / (w + s) Dim y: y = e(l + 1) Dim z: z = 0.0 Dim t, u For k = l To h - 1 If Abs(x) >= Abs(y) Then t = -y / x u = 1 / Sqr(t * t + 1.0) s = t * u Else t = -x / y s = 1 / Sqr(t * t + 1.0) If t < 0 Then s = -s End If u = t * s End If w = d(k) - d(k + 1) t = (w * s + 2 * u * e(k + 1)) * s d(k ) = d(k ) - t d(k + 1) = d(k + 1) + t e(k ) = u * e(k) - s * z e(k + 1) = e(k + 1) * (u * u - s * s) + w * s * u '次の行は固有ベクトルを求めないなら不要 For i = 0 To N x = a(k )(i) y = a(k + 1)(i) a(k )(i) = u * x - s * y a(k + 1)(i) = s * x + u * y Next If k < N - 1 Then x = e(k + 1) y = -s * e(k + 2) z = y e(k + 2) = u * e(k + 2) End If Next WScript.StdOut.Write Right(Space(3) & j, 3) & vbTab Call disp_vector(d) '収束判定 If (Abs(e(h)) < 0.00000000001) Then Exit For Next e(0) = 1.0 Do While (Abs(e(h)) < 0.00000000001) h = h - 1 Loop Loop '次の行は固有ベクトルを求めないなら不要 For k = 0 To N - 1 l = k For i = k + 1 To N If d(i) > d(l) Then l = i End If Next t = d(k) d(k) = d(l) d(l) = t For i = 0 To N t = a(k)(i) a(k)(i) = a(l)(i) a(l)(i) = t Next Next End Sub '１次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '２次元配列を表示 Private Sub disp_matrix(ByVal matrix()) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub
Z:\>cscript //nologo Z:\1106.vbs 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.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 d = [1.0 ,0.0 ,0.0 ,0.0] var e = [1.0 ,0.0 ,0.0 ,0.0] // ハウスホルダー変換 tridiagonalize(a, d, e) // QR分解 decomp(a, d, e) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") disp_vector(d) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvector") disp_matrix(a) } // ハウスホルダー変換 function tridiagonalize(a, d, e) { var v = [0.0 ,0.0 ,0.0 ,0.0] for (var k = 0; k < N - 2; k++) { var w = [0.0 ,0.0 ,0.0 ,0.0] d[k] = a[k][k] var t = 0.0 for (var i = k + 1; i < N; i++) { w[i] = a[i][k] t += w[i] * w[i] } var s = Math.sqrt(t) if (w[k + 1] < 0) s = -s if (Math.abs(s) < 0.00000000001) e[k + 1] = 0.0 else { e[k + 1] = -s w[k + 1] += s s = 1 / Math.sqrt(w[k + 1] * s) for (var i = k + 1; i < N; i++) w[i] *= s for (var i = k + 1; i < N; i++) { s = 0.0 for (var j = k + 1; j < N; j++) { if (j <= i) s += a[i][j] * w[j] else s += a[j][i] * w[j] } v[i] = s } s = 0.0 for (var i = k + 1; i < N; i++) s += w[i] * v[i] s /= 2.0 for (var i = k + 1; i < N; i++) v[i] -= s * w[i] for (var i = k + 1; i < N; i++) for (var j = k + 1; j <= i; j++) a[i][j] -= w[i] * v[j] + w[j] * v[i] // 次の行は固有ベクトルを求めないなら不要 for (var i = k + 1; i < N; i++) a[i][k] = w[i] } } d[N - 2] = a[N - 2][N - 2] d[N - 1] = a[N - 1][N - 1] e[0] = 0.0 e[N - 1] = a[N - 1][N - 2] // 次の行は固有ベクトルを求めないなら不要 for (var k = N - 1; k >= 0; k--) { var w = [0.0 ,0.0 ,0.0 ,0.0] if (k < N - 2) { for (var i = k + 1; i < N; i++) w[i] = a[i][k] for (var i = k + 1; i < N; i++) { var s = 0.0 for (var j = k + 1; j < N; j++) s += a[i][j] * w[j] v[i] = s } for (var i = k + 1; i < N; i++) for (var j = k + 1; j < N; j++) a[i][j] -= v[i] * w[j] } for (var i = 0; i < N; i++) a[i][k] = 0.0 a[k][k] = 1.0 } } // QR分解 function decomp(a, d, e) { e[0] = 1.0 var h = N - 1 while (Math.abs(e[h]) < 0.00000000001) h-- while (h > 0) { e[0] = 0.0 var l = h - 1 while (Math.abs(e[l]) >= 0.00000000001) l-- for (var j = 1; j <= 100; j++) { var w = (d[h - 1] - d[h]) / 2.0 var s = Math.sqrt(w * w + e[h] * e[h]) if (w < 0.0) s = -s var x = d[l] - d[h] + e[h] * e[h] / (w + s) var y = e[l + 1] var z = 0.0 var t = 0.0 var u = 0.0 for (var k = l; k < h; k++) { if (Math.abs(x) >= Math.abs(y)) { t = -y / x u = 1 / Math.sqrt(t * t + 1.0) s = t * u } else { t = -x / y s = 1 / Math.sqrt(t * t + 1.0) if (t < 0) s = -s u = t * s } w = d[k] - d[k + 1] t = (w * s + 2 * u * e[k + 1]) * s d[k ] = d[k ] - t d[k + 1] = d[k + 1] + t e[k ] = u * e[k] - s * z e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for (var i = 0; i < N; i++) { x = a[k ][i] y = a[k + 1][i] a[k ][i] = u * x - s * y a[k + 1][i] = s * x + u * y } if (k < N - 2) { x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] } } WScript.StdOut.Write((" " + j).slice( -3) + "\t") disp_vector(d) // 収束判定 if (Math.abs(e[h]) < 0.00000000001) break } e[0] = 1.0 while (Math.abs(e[h]) < 0.00000000001) h-- } // 次の行は固有ベクトルを求めないなら不要 for (var k = 0; k < N - 1; k++) { var l = k for (var i = k + 1; i < N; i++) if (d[i] > d[l]) l = i var t = d[k] d[k] = d[l] d[l] = t for (var i = 0; i < N; i++) { t = a[k][i] a[k][i] = a[l][i] a[l][i] = t } } } // １次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // ２次元配列を表示 function disp_matrix(matrix) { for (var row = 0; row < N; row++) { for (var col = 0; col < N; col++) WScript.StdOut.Write((" " + matrix[row][col].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }
Z:\>cscript //nologo Z:\1106.js 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
PowerShell
set-variable -name N -value 3 -option constant # １次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # ２次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } # ハウスホルダー変換 function tridiagonalize($a, $d, $e) { $v = (0.0 ,0.0 ,0.0 ,0.0) foreach ($k in 0..($N - 2)) { $d[$k] = $a[$k][$k] $t = 0.0 $w = (0.0, 0.0, 0.0, 0.0) foreach ($i in ($k + 1)..$N) { $w[$i] = $a[$i][$k] $t += ($w[$i] * $w[$i]) } $s = [Math]::Sqrt($t) if ($w[$k + 1] -lt 0) { $s = -$s } if ([Math]::Abs($s) -lt 0.00000000001) { $e[$k + 1] = 0.0 } else { $e[$k + 1] = -$s $w[$k + 1] += $s $s = 1 / [Math]::Sqrt($w[$k + 1] * $s) foreach ($i in ($k + 1)..$N) { $w[$i] *= $s } foreach ($i in ($k + 1)..$N) { $s = 0.0 foreach ($j in ($k + 1)..$N) { if ($j -le $i) { $s += $a[$i][$j] * $w[$j] } else { $s += $a[$j][$i] * $w[$j] } } $v[$i] = $s } $s = 0 foreach ($i in ($k + 1)..$N) { $s += $w[$i] * $v[$i] } $s /= 2 foreach ($i in ($k + 1)..$N) { $v[$i] -= $s * $w[$i] } foreach ($i in ($k + 1)..$N) { foreach ($j in ($k + 1)..$i) { $a[$i][$j] -= $w[$i] * $v[$j] + $w[$j] * $v[$i] } } # 次の行は固有ベクトルを求めないなら不要 foreach ($i in ($k + 1)..$N) { $a[$i][$k] = $w[$i] } } } $d[$N - 1] = $a[$N - 1][$N - 1] $d[$N] = $a[$N][$N] $e[0] = 0.0 $e[$N] = $a[$N][$N - 1] # 次の行は固有ベクトルを求めないなら不要 foreach ($k in ($N..0)) { $w = (0.0, 0.0, 0.0, 0.0) if ($k -lt ($N - 1)) { foreach ($i in ($k + 1)..$N) { $w[$i] = $a[$i][$k] } foreach ($i in ($k + 1)..$N) { $s = 0.0 foreach ($j in ($k + 1)..$N) { $s += $a[$i][$j] * $w[$j] } $v[$i] = $s } foreach ($i in ($k + 1)..$N) { foreach ($j in ($k + 1)..$N) { $a[$i][$j] -= $v[$i] * $w[$j] } } } foreach ($i in 0..$N) { $a[$i][$k] = 0.0 } $a[$k][$k] = 1.0 } } # QR分解 function decomp($a, $d, $e) { $e[0] = 1.0 $h = $N while ([Math]::Abs($e[$h]) -lt 0.00000000001) { $h-- } while ($h -gt 0) { $e[0] = 0.0 $l = $h - 1 while ([Math]::Abs($e[$l]) -ge 0.00000000001) { $l-- } foreach ($j in 1..100) { $w = ($d[$h - 1] - $d[$h]) / 2.0 $s = [Math]::Sqrt($w * $w + $e[$h] * $e[$h]) if ($w -lt 0.0) { $s = - $s } $x = $d[$l] - $d[$h] + $e[$h] * $e[$h] / ($w + $s) $y = $e[$l + 1] $z = 0.0 foreach ($k in $l..($h - 1)) { if ([Math]::Abs($x) -ge [Math]::Abs($y)) { $t = -$y / $x $u = 1 / [Math]::Sqrt($t * $t + 1.0) $s = $t * $u } else { $t = -$x / $y $s = 1 / [Math]::Sqrt($t * $t + 1.0) if ($t -lt 0) { $s = -$s } $u = $t * $s } $w = $d[$k] - $d[$k + 1] $t = ($w * $s + 2 * $u * $e[$k + 1]) * $s $d[$k ] = $d[$k ] - $t $d[$k + 1] = $d[$k + 1] + $t $e[$k ] = $u * $e[$k] - $s * $z $e[$k + 1] = $e[$k + 1] * ($u * $u - $s * $s) + $w * $s * $u # 次の行は固有ベクトルを求めないなら不要 foreach ($i in 0..$N) { $x = $a[$k ][$i] $y = $a[$k + 1][$i] $a[$k ][$i] = $u * $x - $s * $y $a[$k + 1][$i] = $s * $x + $u * $y } if ($k -lt $N - 1) { $x = $e[$k + 1] $y = -$s * $e[$k + 2] $z = $y $e[$k + 2] = $u * $e[$k + 2] } } Write-Host ([String]::Format("{0,3:D}`t", $j)) -nonewline disp_vector $d # 収束判定 if ([Math]::Abs($e[$h]) -lt 0.00000000001) { break } } $e[0] = 1.0 while ([Math]::Abs($e[$h]) -lt 0.00000000001) { $h-- } } # 次の行は固有ベクトルを求めないなら不要 foreach ($k in 0..($N - 1)) { $l = $k foreach ($i in ($k + 1)..$N) { if ($d[$i] -gt $d[$l]) { $l = $i } } $t = $d[$k] $d[$k] = $d[$l] $d[$l] = $t foreach ($i in 0..$N) { $t = $a[$k][$i] $a[$k][$i] = $a[$l][$i] $a[$l][$i] = $t } } } # ハウスホルダー変換と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) $d = (0.0 ,0.0 ,0.0 ,0.0) $e = (0.0 ,0.0 ,0.0 ,0.0) # ハウスホルダー変換 tridiagonalize $a $d $e # QR分解 decomp $a $d $e Write-Host Write-Host "固有値" disp_vector $d Write-Host Write-Host "固有ベクトル" disp_matrix $a
Z:\>powershell -file Z:\1106.ps1 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 固有値 10.0000000000 5.0000000000 2.0000000000 1.0000000000 固有ベクトル 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.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 @d = (0.0, 0.0, 0.0, 0.0); my @e = (0.0, 0.0, 0.0, 0.0); # ハウスホルダー変換 tridiagonalize(\@a, \@d, \@e); # QR分解 decomp(\@a, \@d, \@e); print "\neigenvalue\n"; disp_vector(\@d); print "\neigenvector\n"; disp_matrix(\@a); } # ハウスホルダー変換 sub tridiagonalize { my ($a, $d, $e) = @_; my @v = (0.0, 0.0, 0.0, 0.0); for $k (0..(N - 2)) { my @w = (0.0, 0.0, 0.0, 0.0); $$d[$k] = $$a[$k][$k]; my $t = 0.0; for $i (($k + 1)..N) { $w[$i] = $$a[$i][$k]; $t += ($w[$i] * $w[$i]); } my $s = sqrt($t); $s = -$s if ($w[$k + 1] < 0); if (abs($s) < 0.00000000001) { $$e[$k + 1] = 0.0; } else { $$e[$k + 1] = -$s; $w[$k + 1] += $s; $s = 1 / sqrt($w[$k + 1] * $s); for $i (($k + 1)..N) { $w[$i] *= $s; } for $i (($k + 1)..N) { $s = 0.0; for $j (($k + 1)..N) { if ($j <= $i) { $s += $$a[$i][$j] * $w[$j]; } else { $s += $$a[$j][$i] * $w[$j]; } } $v[$i] = $s; } $s = 0; for $i (($k + 1)..N) { $s += $w[$i] * $v[$i]; } $s /= 2; for $i (($k + 1)..N) { $v[$i] -= $s * $w[$i]; } for $i (($k + 1)..N) { for $j (($k + 1)..$i) { $$a[$i][$j] -= $w[$i] * $v[$j] + $w[$j] * $v[$i]; } } # 次の行は固有ベクトルを求めないなら不要 for $i (($k + 1)..N) { $$a[$i][$k] = $w[$i]; } } } $$d[N - 1] = $$a[N - 1][N - 1]; $$d[N] = $$a[N][N]; $$e[0] = 0.0; $$e[N] = $$a[N][N - 1]; # 次の行は固有ベクトルを求めないなら不要 for $k (reverse(0..N)) { $w = (0.0, 0.0, 0.0, 0.0); if ($k < (N - 1)) { for $i (($k + 1)..N) { $w[$i] = $$a[$i][$k]; } for $i (($k + 1)..N) { $s = 0.0; for $j (($k + 1)..N) { $s += $$a[$i][$j] * $w[$j]; } $v[$i] = $s; } for $i (($k + 1)..N) { for $j (($k + 1)..N) { $$a[$i][$j] -= $v[$i] * $w[$j]; } } } for $i (0..N) { $$a[$i][$k] = 0.0; } $$a[$k][$k] = 1.0; } } # QR分解 sub decomp { my ($a, $d, $e) = @_; $$e[0] = 1.0; my $h = N; $h-- while (abs($$e[$h]) < 0.00000000001); while ($h > 0) { $$e[0] = 0.0; my $l = $h - 1; $l-- while (abs($$e[$l]) >= 0.00000000001); for $j (1..100) { my $w = ($$d[$h - 1] - $$d[$h]) / 2.0; my $s = sqrt($w * $w + $$e[$h] * $$e[$h]); $s = - $s if ($w < 0.0); my $x = $$d[$l] - $$d[$h] + $$e[$h] * $$e[$h] / ($w + $s); my $y = $$e[$l + 1]; my $z = 0.0; my $t = 0.0; my $u = 0.0; for $k ($l..($h - 1)) { if (abs($x) >= abs($y)) { $t = -$y / $x; $u = 1 / sqrt($t * $t + 1.0); $s = $t * $u; } else { $t = -$x / $y; $s = 1 / sqrt($t * $t + 1.0); $s = -$s if ($t < 0.0); $u = $t * $s; } $w = $$d[$k] - $$d[$k + 1]; $t = ($w * $s + 2 * $u * $$e[$k + 1]) * $s; $$d[$k ] = $$d[$k ] - $t; $$d[$k + 1] = $$d[$k + 1] + $t; $$e[$k ] = $u * $$e[$k] - $s * $z; $$e[$k + 1] = $$e[$k + 1] * ($u * $u - $s * $s) + $w * $s * $u; # 次の行は固有ベクトルを求めないなら不要 for $i (0..N) { $x = $$a[$k ][$i]; $y = $$a[$k + 1][$i]; $$a[$k ][$i] = $u * $x - $s * $y; $$a[$k + 1][$i] = $s * $x + $u * $y; } if ($k < N - 1) { $x = $$e[$k + 1]; $y = -$s * $$e[$k + 2]; $z = $y; $$e[$k + 2] = $u * $$e[$k + 2]; } } printf("%3d\t", $j); disp_vector(\@$d); # 収束判定 last if (abs($$e[$h]) < 0.00000000001); } $$e[0] = 1.0; $h-- while (abs($$e[$h]) < 0.00000000001); } # 次の行は固有ベクトルを求めないなら不要 for $k (0..(N - 1)) { $l = $k; for $i (($k + 1)..N) { $l = $i if ($$d[$i] > $$d[$l]); } $t = $$d[$k]; $$d[$k] = $$d[$l]; $$d[$l] = $t; for $i (0..N) { $t = $$a[$k][$i]; $$a[$k][$i] = $$a[$l][$i]; $$a[$l][$i] = $t; } } } # １次元配列を表示 sub disp_vector { my ($x) = @_; for $i (0..N) { printf("%14.10f\t", $$x[$i]); } print "\n"; } # ２次元配列を表示 sub disp_matrix { my ($matrix) = @_; foreach $row (@$matrix) { foreach $col (@$row) { printf("%14.10f\t", $col); } print "\n"; } }
Z:\>perl Z:\1106.pl 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.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]]; $d = [0.0, 0.0, 0.0, 0.0]; $e = [0.0, 0.0, 0.0, 0.0]; # ハウスホルダー変換 tridiagonalize($a, $d, $e); # QR分解 decomp($a, $d, $e); print "\neigenvalue\n"; disp_vector($d); print "\neigenvector\n"; disp_matrix($a); } # ハウスホルダー変換 function tridiagonalize(&$a, &$d, &$e) { $v = [0.0, 0.0, 0.0, 0.0]; foreach (range(0, (N - 2)) as $k) { $w = [0.0, 0.0, 0.0, 0.0]; $d[$k] = $a[$k][$k]; $t = 0.0; foreach (range(($k + 1), N) as $i) { $w[$i] = $a[$i][$k]; $t += ($w[$i] * $w[$i]); } $s = sqrt($t); if ($w[$k + 1] < 0) $s = -$s; if (abs($s) < 0.00000000001) { $e[$k + 1] = 0.0; } else { $e[$k + 1] = -$s; $w[$k + 1] += $s; $s = 1 / sqrt($w[$k + 1] * $s); foreach (range(($k + 1), N) as $i) { $w[$i] *= $s; } foreach (range(($k + 1), N) as $i) { $s = 0.0; foreach (range(($k + 1), N) as $j) { if ($j <= $i) { $s += $a[$i][$j] * $w[$j]; } else { $s += $a[$j][$i] * $w[$j]; } } $v[$i] = $s; } $s = 0; foreach (range(($k + 1), N) as $i) { $s += $w[$i] * $v[$i]; } $s /= 2; foreach (range(($k + 1), N) as $i) { $v[$i] -= $s * $w[$i]; } foreach (range(($k + 1), N) as $i) { foreach (range(($k + 1), $i) as $j) { $a[$i][$j] -= $w[$i] * $v[$j] + $w[$j] * $v[$i]; } } # 次の行は固有ベクトルを求めないなら不要 foreach (range(($k + 1), N) as $i) { $a[$i][$k] = $w[$i]; } } } $d[N - 1] = $a[N - 1][N - 1]; $d[N] = $a[N][N]; $e[0] = 0.0; $e[N] = $a[N][N - 1]; # 次の行は固有ベクトルを求めないなら不要 foreach (range(N, 0) as $k) { $w = [0.0, 0.0, 0.0, 0.0]; if ($k < (N - 1)) { foreach (range(($k + 1), N) as $i) { $w[$i] = $a[$i][$k]; } foreach (range(($k + 1), N) as $i) { $s = 0.0; foreach (range(($k + 1), N) as $j) { $s += $a[$i][$j] * $w[$j]; } $v[$i] = $s; } foreach (range(($k + 1), N) as $i) { foreach (range(($k + 1), N) as $j) { $a[$i][$j] -= $v[$i] * $w[$j]; } } } foreach (range(0, N) as $i) { $a[$i][$k] = 0.0; } $a[$k][$k] = 1.0; } } # QR分解 function decomp(&$a, &$d, &$e) { $e[0] = 1.0; $h = N; while (abs($e[$h]) < 0.00000000001) $h--; while ($h > 0) { $e[0] = 0.0; $l = $h - 1; while (abs($e[$l]) >= 0.00000000001) $l--; foreach (range(1, 100) as $j) { $w = ($d[$h - 1] - $d[$h]) / 2.0; $s = sqrt($w * $w + $e[$h] * $e[$h]); if ($w < 0.0) $s = - $s; $x = $d[$l] - $d[$h] + $e[$h] * $e[$h] / ($w + $s); $y = $e[$l + 1]; $z = 0.0; $t = 0.0; $u = 0.0; foreach (range(l, ($h - 1)) as $k) { if (abs($x) >= abs($y)) { $t = -$y / $x; $u = 1 / sqrt($t * $t + 1.0); $s = $t * $u; } else { $t = -$x / $y; $s = 1 / sqrt($t * $t + 1.0); if ($t < 0.0) $s = -$s; $u = $t * $s; } $w = $d[$k] - $d[$k + 1]; $t = ($w * $s + 2 * $u * $e[$k + 1]) * $s; $d[$k ] = $d[$k ] - $t; $d[$k + 1] = $d[$k + 1] + $t; $e[$k ] = $u * $e[$k] - $s * $z; $e[$k + 1] = $e[$k + 1] * ($u * $u - $s * $s) + $w * $s * $u; # 次の行は固有ベクトルを求めないなら不要 foreach (range(0, N) as $i) { $x = $a[$k ][$i]; $y = $a[$k + 1][$i]; $a[$k ][$i] = $u * $x - $s * $y; $a[$k + 1][$i] = $s * $x + $u * $y; } if ($k < N - 1) { $x = $e[$k + 1]; $y = -$s * $e[$k + 2]; $z = $y; $e[$k + 2] = $u * $e[$k + 2]; } } printf("%3d\t", $j); disp_vector($d); # 収束判定 if (abs($e[$h]) < 0.00000000001) break; } $e[0] = 1.0; while (abs($e[$h]) < 0.00000000001) $h--; } # 次の行は固有ベクトルを求めないなら不要 foreach (range(0, (N - 1)) as $k) { $l = $k; foreach (range(($k + 1), N) as $i) { if ($d[$i] > $d[$l]) $l = $i; } $t = $d[$k]; $d[$k] = $d[$l]; $d[$l] = $t; foreach (range(0, N) as $i) { $t = $a[$k][$i]; $a[$k][$i] = $a[$l][$i]; $a[$l][$i] = $t; } } } # １次元配列を表示 function disp_vector($row) { foreach ($row as $col) { printf("%14.10f\t", $col); } print "\n"; } # ２次元配列を表示 function disp_matrix($matrix) { foreach ($matrix as $row) { foreach ($row as $col) { printf("%14.10f\t", $col); } print "\n"; } } ?>
Z:\>php Z:\1106.php 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Python
# coding: Shift_JIS import math N = 4 # １次元配列を表示 def disp_vector(row): for col in row: print "%14.10f\t" % col, print "" # ２次元配列を表示 def disp_matrix(matrix): for row in matrix: for col in row: print "%14.10f\t" % col, print "" # ハウスホルダー変換 def tridiagonalize(a, d, e): v = [0.0 ,0.0 ,0.0 ,0.0] for k in range(0, N - 2, 1): w = [0.0 ,0.0 ,0.0 ,0.0] d[k] = a[k][k] t = 0.0 for i in range(k + 1, N, 1): w[i] = a[i][k] t += w[i] * w[i] s = math.sqrt(t) if (w[k + 1] < 0): s = -s if (abs(s) < 0.00000000001): e[k + 1] = 0.0 else: e[k + 1] = -s w[k + 1] += s s = 1 / math.sqrt(w[k + 1] * s) for i in range(k + 1, N, 1): w[i] *= s for i in range(k + 1, N, 1): s = 0.0 for j in range(k + 1, N, 1): if (j <= i): s += a[i][j] * w[j] else: s += a[j][i] * w[j] v[i] = s s = 0.0 for i in range(k + 1, N, 1): s += w[i] * v[i] s /= 2.0 for i in range(k + 1, N, 1): v[i] -= s * w[i] for i in range(k + 1, N, 1): for j in range(k + 1, i + 1, 1): a[i][j] -= w[i] * v[j] + w[j] * v[i] # 次の行は固有ベクトルを求めないなら不要 for i in range(k + 1, N, 1): a[i][k] = w[i] d[N - 2] = a[N - 2][N - 2] d[N - 1] = a[N - 1][N - 1] e[0] = 0.0 e[N - 1] = a[N - 1][N - 2] # 次の行は固有ベクトルを求めないなら不要 for k in range(N - 1, -1, -1): w = [0.0 ,0.0 ,0.0 ,0.0] if (k < N - 2): for i in range(k + 1, N, 1): w[i] = a[i][k] for i in range(k + 1, N, 1): s = 0.0 for j in range(k + 1, N, 1): s += a[i][j] * w[j] v[i] = s for i in range(k + 1, N, 1): for j in range(k + 1, N, 1): a[i][j] -= v[i] * w[j] for i in range(0, N, 1): a[i][k] = 0.0 a[k][k] = 1.0 # QR分解 def decomp(a, d, e): e[0] = 1.0 h = N - 1 while (abs(e[h]) < 0.00000000001): h -= 1 while (h > 0): e[0] = 0.0 l = h - 1 while (abs(e[l]) >= 0.00000000001): l -= 1 for j in range(1, 100, 1): w = (d[h - 1] - d[h]) / 2.0 s = math.sqrt(w * w + e[h] * e[h]) if (w < 0.0): s = -s x = d[l] - d[h] + e[h] * e[h] / (w + s) y = e[l + 1] z = 0.0 t = 0.0 u = 0.0 for k in range(l, h, 1): if (abs(x) >= abs(y)): t = -y / x u = 1 / math.sqrt(t * t + 1.0) s = t * u else: t = -x / y s = 1 / math.sqrt(t * t + 1.0) if (t < 0): s = -s u = t * s w = d[k] - d[k + 1] t = (w * s + 2 * u * e[k + 1]) * s d[k ] = d[k ] - t d[k + 1] = d[k + 1] + t e[k ] = u * e[k] - s * z e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u # 次の行は固有ベクトルを求めないなら不要 for i in range(0, N, 1): x = a[k ][i] y = a[k + 1][i] a[k ][i] = u * x - s * y a[k + 1][i] = s * x + u * y if (k < N - 2): x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] print "%3d\t" % j, disp_vector(d) # 収束判定 if (abs(e[h]) < 0.00000000001): break e[0] = 1.0 while (abs(e[h]) < 0.00000000001): h -= 1 # 次の行は固有ベクトルを求めないなら不要 for k in range(0, N - 1, 1): l = k for i in range(k + 1, N, 1): if (d[i] > d[l]): l = i t = d[k] d[k] = d[l] d[l] = t for i in range(0, N, 1): t = a[k][i] a[k][i] = a[l][i] a[l][i] = t # ハウスホルダー変換と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]] d = [1.0 ,0.0 ,0.0 ,0.0] e = [1.0 ,0.0 ,0.0 ,0.0] # ハウスホルダー変換 tridiagonalize(a, d, e) # QR分解 decomp(a, d, e) print "" print "eigenvalue" disp_vector(d) print "" print "eigenvector" disp_matrix(a)
Z:\>python Z:\1106.py 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Ruby
# coding: Shift_JIS N = 3 # １次元配列を表示 def disp_vector(row) row.each do |col| printf("%14.10f\t", col) end puts "" end # ２次元配列を表示 def disp_matrix(matrix) matrix.each do |row| row.each do |col| printf("%14.10f\t", col) end puts "" end end # ハウスホルダー変換 def tridiagonalize(a, d, e) v = [0.0 ,0.0 ,0.0 ,0.0] (0..(N - 2)).each do |k| w = [0.0 ,0.0 ,0.0 ,0.0] d[k] = a[k][k] t = 0.0 ((k + 1)..N).each do |i| w[i] = a[i][k] t += w[i] * w[i] end s = Math.sqrt(t) s = -s if (w[k + 1] < 0) if (s.abs < 0.00000000001) e[k + 1] = 0.0 else e[k + 1] = -s w[k + 1] += s s = 1 / Math.sqrt(w[k + 1] * s) ((k + 1)..N).each do |i| w[i] *= s end ((k + 1)..N).each do |i| s = 0.0 ((k + 1)..N).each do |j| if (j <= i) s += a[i][j] * w[j] else s += a[j][i] * w[j] end end v[i] = s end s = 0.0 ((k + 1)..N).each do |i| s += w[i] * v[i] end s /= 2.0 ((k + 1)..N).each do |i| v[i] -= s * w[i] end ((k + 1)..N).each do |i| ((k + 1)..i).each do |j| a[i][j] -= w[i] * v[j] + w[j] * v[i] end end # 次の行は固有ベクトルを求めないなら不要 ((k + 1)..N).each do |i| a[i][k] = w[i] end end end d[N - 1] = a[N - 1][N - 1] d[N] = a[N][N] e[0] = 0.0 e[N] = a[N][N - 1] # 次の行は固有ベクトルを求めないなら不要 (0..N).reverse_each do |k| w = [0.0 ,0.0 ,0.0 ,0.0] if (k < N - 1) ((k + 1)..N).each do |i| w[i] = a[i][k] end ((k + 1)..N).each do |i| s = 0.0 ((k + 1)..N).each do |j| s += a[i][j] * w[j] end v[i] = s end ((k + 1)..N).each do |i| ((k + 1)..N).each do |j| a[i][j] -= v[i] * w[j] end end end (0..N).each do |i| a[i][k] = 0.0 end a[k][k] = 1.0 end end # QR分解 def decomp(a, d, e) e[0] = 1.0 h = N - 1 h -= 1 while ((e[h]).abs < 0.00000000001) while (h > 0) e[0] = 0.0 l = h - 1 l -= 1 while ((e[l]).abs >= 0.00000000001) (1..100).each do |j| w = (d[h - 1] - d[h]) / 2.0 s = Math.sqrt(w * w + e[h] * e[h]) s = -s if (w < 0.0) x = d[l] - d[h] + e[h] * e[h] / (w + s) y = e[l + 1] z = 0.0 t = 0.0 u = 0.0 (l..(h - 1)).each do |k| if (x.abs >= y.abs) t = -y / x u = 1 / Math.sqrt(t * t + 1.0) s = t * u else t = -x / y s = 1 / Math.sqrt(t * t + 1.0) s = -s if (t < 0) u = t * s end w = d[k] - d[k + 1] t = (w * s + 2 * u * e[k + 1]) * s d[k ] = d[k ] - t d[k + 1] = d[k + 1] + t e[k ] = u * e[k] - s * z e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u # 次の行は固有ベクトルを求めないなら不要 (0..N).each do |i| x = a[k ][i] y = a[k + 1][i] a[k ][i] = u * x - s * y a[k + 1][i] = s * x + u * y end if (k < N - 1) x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] end end printf("%3d\t", j) disp_vector(d) # 収束判定 break if ((e[h]).abs < 0.00000000001) end e[0] = 1.0 h -= 1 while ((e[h]).abs < 0.00000000001) end # 次の行は固有ベクトルを求めないなら不要 (0..(N - 1)).each do |k| l = k ((k + 1)..N).each do |i| l = i if (d[i] > d[l]) end t = d[k] d[k] = d[l] d[l] = t (0..N).each do |i| t = a[k][i] a[k][i] = a[l][i] a[l][i] = t 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]] d = [1.0 ,0.0 ,0.0 ,0.0] e = [1.0 ,0.0 ,0.0 ,0.0] # ハウスホルダー変換 tridiagonalize(a, d, e) # QR分解 decomp(a, d, e) puts "" puts "eigenvalue" disp_vector(d) puts "" puts "eigenvector" disp_matrix(a)
Z:\>ruby Z:\1106.rb 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.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 d = [1.0 ,0.0 ,0.0 ,0.0] as double[] def e = [1.0 ,0.0 ,0.0 ,0.0] as double[] // ハウスホルダー変換 tridiagonalize(a, d, e) // QR分解 decomp(a, d, e) println() println("eigenvalue") disp_vector(d) println() println("eigenvector") disp_matrix(a) } // ハウスホルダー変換 def tridiagonalize(a, d, e) { def v = [0.0 ,0.0 ,0.0 ,0.0] as double[] for (k in 0..(N - 2)) { def w = [0.0 ,0.0 ,0.0 ,0.0] as double[] d[k] = a[k][k] t = 0.0 for (i in (k + 1)..N) { w[i] = a[i][k] t += w[i] * w[i] } s = Math.sqrt(t) if (w[k + 1] < 0) s = -s if (Math.abs(s) < 0.00000000001) e[k + 1] = 0.0 else { e[k + 1] = -s w[k + 1] += s s = 1 / Math.sqrt(w[k + 1] * s) for (i in (k + 1)..N) w[i] *= s for (i in (k + 1)..N) { s = 0.0 for (j in (k + 1)..N) { if (j <= i) s += a[i][j] * w[j] else s += a[j][i] * w[j] } v[i] = s } s = 0.0 for (i in (k + 1)..N) s += w[i] * v[i] s /= 2.0 for (i in (k + 1)..N) v[i] -= s * w[i] for (i in (k + 1)..N) for (j in (k + 1)..i) a[i][j] -= w[i] * v[j] + w[j] * v[i] // 次の行は固有ベクトルを求めないなら不要 for (i in (k + 1)..N) a[i][k] = w[i] } } d[N - 1] = a[N - 1][N - 1] d[N] = a[N][N] e[0] = 0.0 e[N] = a[N][N - 1] // 次の行は固有ベクトルを求めないなら不要 for (k in N..0) { w = [0.0 ,0.0 ,0.0 ,0.0] if (k < N - 1) { for (i in (k + 1)..N) w[i] = a[i][k] for (i in (k + 1)..N) { s = 0.0 for (j in (k + 1)..N) s += a[i][j] * w[j] v[i] = s } for (i in (k + 1)..N) for (j in (k + 1)..N) a[i][j] -= v[i] * w[j] } for (i in 0..N) a[i][k] = 0.0 a[k][k] = 1.0 } } // QR分解 def decomp(a, d, e) { e[0] = 1.0 h = N - 1 while (Math.abs(e[h]) < 0.00000000001) h-- while (h > 0) { e[0] = 0.0 l = h - 1 while (Math.abs(e[l]) >= 0.00000000001) l-- for (j in 1..100) { w = (d[h - 1] - d[h]) / 2.0 s = Math.sqrt(w * w + e[h] * e[h]) if (w < 0.0) s = -s x = d[l] - d[h] + e[h] * e[h] / (w + s) y = e[l + 1] z = 0.0 t = 0.0 u = 0.0 for (k in l..(h - 1)) { if (Math.abs(x) >= Math.abs(y)) { t = -y / x u = 1 / Math.sqrt(t * t + 1.0) s = t * u } else { t = -x / y s = 1 / Math.sqrt(t * t + 1.0) if (t < 0) s = -s u = t * s } w = d[k] - d[k + 1] t = (w * s + 2 * u * e[k + 1]) * s d[k ] = d[k ] - t d[k + 1] = d[k + 1] + t e[k ] = u * e[k] - s * z e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for (i in 0..N) { x = a[k ][i] y = a[k + 1][i] a[k ][i] = u * x - s * y a[k + 1][i] = s * x + u * y } if (k < N - 1) { x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] } } printf("%3d\t" , j) disp_vector(d) // 収束判定 if (Math.abs(e[h]) < 0.00000000001) break } e[0] = 1.0 while (Math.abs(e[h]) < 0.00000000001) h-- } // 次の行は固有ベクトルを求めないなら不要 for (k in 0..(N - 1)) { l = k for (i in (k + 1)..N) if (d[i] > d[l]) l = i t = d[k] d[k] = d[l] d[l] = t for (i in 0..N) { t = a[k][i] a[k][i] = a[l][i] a[l][i] = t } } } // １次元配列を表示 def disp_vector(row) { for (col in row) { printf("%14.10f\t" , col) } println() } // ２次元配列を表示 def disp_matrix(matrix) { for (row in matrix) { for (col in row) { printf("%14.10f\t" , col) } println() } }
Z:\>gvy Gvy1106.gvy 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Pascal
program Pas1106(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // １次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // ２次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var row, col:Integer; begin for row := 0 to N do begin for col := 0 to N do write(format('%14.10f'#9, [matrix[row][col]])); writeln(); end; end; // ハウスホルダー変換 procedure tridiagonalize(var a:TwoDimArray; var d:array of Double; var e:array of Double); var w:array [0..N] of Double; v:array [0..N] of Double; i, j, k:Integer; s, t:Double; begin for k := 0 to N do v[k] := 0.0; for k := 0 to (N - 2) do begin for i := 0 to N do w[i] := 0.0; d[k] := a[k][k]; t := 0.0; for i := (k + 1) to N do begin w[i] := a[i][k]; t := t + w[i] * w[i]; end; s := Sqrt(t); if w[k + 1] < 0 then s := -s; if Abs(s) < 0.00000000001 then begin e[k + 1] := 0.0; end else begin e[k + 1] := -s; w[k + 1] := w[k + 1] + s; s := 1.0 / Sqrt(w[k + 1] * s); for i := (k + 1) to N do w[i] := w[i] * s; for i := (k + 1) to N do begin s := 0.0; for j := (k + 1) to N do if j <= i then s := s + a[i][j] * w[j] else s := s + a[j][i] * w[j]; v[i] := s; end; s := 0.0; for i := (k + 1) to N do s := s + w[i] * v[i]; s := s / 2.0; for i := (k + 1) to N do v[i] := v[i] - s * w[i]; for i := (k + 1) to N do for j := (k + 1) to i do a[i][j] := a[i][j] - (w[i] * v[j] + w[j] * v[i]); // 次の行は固有ベクトルを求めないなら不要 for i := (k + 1) to N do a[i][k] := w[i]; end; end; d[N - 1] := a[N - 1][N - 1]; d[N] := a[N][N]; e[0] := 0.0; e[N] := a[N][N - 1]; // 次の行は固有ベクトルを求めないなら不要 for k := N downto 0 do begin for i := 0 to N do w[i] := 0.0; if k < N - 1 then begin for i := (k + 1) to N do w[i] := a[i][k]; for i := (k + 1) to N do begin s := 0.0; for j := (k + 1) to N do s := s + a[i][j] * w[j]; v[i] := s; end; for i := (k + 1) to N do for j := (k + 1) to N do a[i][j] := a[i][j] - v[i] * w[j]; end; for i := 0 to N do a[i][k] := 0.0; a[k][k] := 1.0; end; end; // QR分解 procedure decomp(var a:TwoDimArray; var d:array of Double; var e:array of Double); var h, l:Integer; i, j, k:Integer; w, s :Double; x, y, z :Double; t, u :Double; begin e[0] := 1.0; h := N; while (Abs(e[h]) < 0.00000000001) do dec(h); while (h > 0) do begin e[0] := 0.0; l := h - 1; while (Abs(e[l]) >= 0.00000000001) do dec(l); for j := 1 to 100 do begin w := (d[h - 1] - d[h]) / 2.0; s := Sqrt(w * w + e[h] * e[h]); if w < 0.0 then s := -s; x := d[l] - d[h] + e[h] * e[h] / (w + s); y := e[l + 1]; z := 0.0; t := 0.0; u := 0.0; for k := l to (h - 1) do begin if Abs(x) >= Abs(y) then begin t := -y / x; u := 1.0 / Sqrt(t * t + 1.0); s := t * u; end else begin t := -x / y; s := 1.0 / Sqrt(t * t + 1.0); if t < 0 then s := -s; u := t * s; end; w := d[k] - d[k + 1]; t := (w * s + 2.0 * u * e[k + 1]) * s; d[k ] := d[k ] - t; d[k + 1] := d[k + 1] + t; e[k ] := u * e[k] - s * z; e[k + 1] := e[k + 1] * (u * u - s * s) + w * s * u; // 次の行は固有ベクトルを求めないなら不要 for i := 0 to N do begin x := a[k , i]; y := a[k + 1, i]; a[k , i] := u * x - s * y; a[k + 1, i] := s * x + u * y; end; if k < N - 1 then begin x := e[k + 1]; y := -s * e[k + 2]; z := y; e[k + 2] := u * e[k + 2]; end; end; write(format('%3d'#9, [j])); disp_vector(d); // 収束判定 if (Abs(e[h]) < 0.00000000001) then break; end; e[0] := 1.0; while (Abs(e[h]) < 0.00000000001) do dec(h); end; // 次の行は固有ベクトルを求めないなら不要 for k := 0 to (N - 1) do begin l := k; for i := (k + 1) to N do if d[i] > d[l] then l := i; t := d[k]; d[k] := d[l]; d[l] := t; for i := 0 to N do begin t := a[k][i]; a[k][i] := a[l][i]; a[l][i] := t; end; 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)); d:array [0..N] of Double; e:array [0..N] of Double; begin // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 decomp(a, d, e); writeln(); writeln('eigenvalue'); disp_vector(d); writeln(); writeln('eigenvector'); disp_matrix(a); end.
Z:\>fpc -v0 -l- Pas1106.pp Z:\>Pas1106 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Ada
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1106 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)); d, e:Long_Float_Array; -- １次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- ２次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for row in 0..N loop for col in 0..N loop Put(matrix(row, col), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- ハウスホルダー変換 procedure tridiagonalize(a:in out Long_Float_TwoDimArray; d:in out Long_Float_Array; e:in out Long_Float_Array) is w, v:Long_Float_Array; t, s:Long_Float; begin for k in 0..N loop v(k) := 0.0; end loop; for k in 0..(N - 2) loop for i in 0..N loop w(i) := 0.0; end loop; d(k) := a(k, k); t := 0.0; for i in (k + 1)..N loop w(i) := a(i, k); t := t + w(i) * w(i); end loop; s := Sqrt(t); if w(k + 1) < 0.0 then s := -s; end if; if Abs(s) < 0.00000000001 then e(k + 1) := 0.0; else e(k + 1) := -s; w(k + 1) := w(k + 1) + s; s := 1.0 / Sqrt(w(k + 1) * s); for i in (k + 1)..N loop w(i) := w(i) * s; end loop; for i in (k + 1)..N loop s := 0.0; for j in (k + 1)..N loop if j <= i then s := s + a(i, j) * w(j); else s := s + a(j, i) * w(j); end if; end loop; v(i) := s; end loop; s := 0.0; for i in (k + 1)..N loop s := s + w(i) * v(i); end loop; s := s / 2.0; for i in (k + 1)..N loop v(i) := v(i) - s * w(i); end loop; for i in (k + 1)..N loop for j in (k + 1)..i loop a(i, j) := a(i, j) - (w(i) * v(j) + w(j) * v(i)); end loop; end loop; -- 次の行は固有ベクトルを求めないなら不要 for i in (k + 1)..N loop a(i, k) := w(i); end loop; end if; end loop; d(N - 1) := a(N - 1, N - 1); d(N) := a(N, N); e(0) := 0.0; e(N) := a(N, N - 1); -- 次の行は固有ベクトルを求めないなら不要 for k in reverse 0..N loop for i in 0..N loop w(i) := 0.0; end loop; if k < N - 1 then for i in (k + 1)..N loop w(i) := a(i, k); end loop; for i in (k + 1)..N loop s := 0.0; for j in (k + 1)..N loop s := s + a(i, j) * w(j); end loop; v(i) := s; end loop; for i in (k + 1)..N loop for j in (k + 1)..N loop a(i, j) := a(i, j) - v(i) * w(j); end loop; end loop; end if; for i in 0..N loop a(i, k) := 0.0; end loop; a(k, k) := 1.0; end loop; end tridiagonalize; -- QR分解 procedure decomp(a:in out Long_Float_TwoDimArray; d:in out Long_Float_Array; e:in out Long_Float_Array) is h, l:Integer; s, t, u, w, x, y, z:Long_Float; begin e(0) := 1.0; h := N; while (Abs(e(h)) < 0.00000000001) loop h := h -1; end loop; while (h > 0) loop e(0) := 0.0; l := h - 1; while (Abs(e(l)) >= 0.00000000001) loop l := l - 1; end loop; for j in 1..100 loop w := (d(h - 1) - d(h)) / 2.0; s := Sqrt(w * w + e(h) * e(h)); if w < 0.0 then s := -s; end if; x := d(l) - d(h) + e(h) * e(h) / (w + s); y := e(l + 1); z := 0.0; t := 0.0; u := 0.0; for k in l..(h - 1) loop if Abs(x) >= Abs(y) then t := -y / x; u := 1.0 / Sqrt(t * t + 1.0); s := t * u; else t := -x / y; s := 1.0 / Sqrt(t * t + 1.0); if t < 0.0 then s := -s; end if; u := t * s; end if; w := d(k) - d(k + 1); t := (w * s + 2.0 * u * e(k + 1)) * s; d(k ) := d(k ) - t; d(k + 1) := d(k + 1) + t; e(k ) := u * e(k) - s * z; e(k + 1) := e(k + 1) * (u * u - s * s) + w * s * u; -- 次の行は固有ベクトルを求めないなら不要 for i in 0..N loop x := a(k , i); y := a(k + 1, i); a(k , i) := u * x - s * y; a(k + 1, i) := s * x + u * y; end loop; if k < N - 1 then x := e(k + 1); y := -s * e(k + 2); z := y; e(k + 2) := u * e(k + 2); end if; end loop; Put(j, Width=> 3); Put(Ascii.HT); disp_vector(d); -- 収束判定 if (Abs(e(h)) < 0.00000000001) then exit; end if; end loop; e(0) := 1.0; while (Abs(e(h)) < 0.00000000001) loop h := h - 1; end loop; end loop; -- 次の行は固有ベクトルを求めないなら不要 for k in 0..(N - 1) loop l := k; for i in (k + 1)..N loop if d(i) > d(l) then l := i; end if; end loop; t := d(k); d(k) := d(l); d(l) := t; for i in 0..N loop t := a(k, i); a(k, i) := a(l, i); a(l, i) := t; end loop; end loop; end decomp; -- ハウスホルダー変換とQR分解で固有値を求める begin -- ハウスホルダー変換; tridiagonalize(a, d, e); -- QR分解 decomp(a, d, e); New_Line; Put_Line("eigenvalue"); disp_vector(d); New_Line; Put_Line("eigenvector"); disp_matrix(a); end Ada1106;
xxxxxx@yyyyyy /Z $ gnatmake Ada1106.adb xxxxxx@yyyyyy /Z $ Ada1106 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 -0.0000000000
VB.NET
Option Explicit Module VB1106 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 d() As Double = New Double(N){} Dim e() As Double = New Double(N){} 'ハウスホルダー変換 tridiagonalize(a, d, e) 'QR分解 decomp(a, d, e) Console.WriteLine() Console.WriteLine("eigenvalue") disp_vector(d) Console.WriteLine() Console.WriteLine("eigenvector") disp_matrix(a) End Sub 'ハウスホルダー変換 Private Sub tridiagonalize(ByRef a(,) As Double, ByRef d() As Double, ByRef e() As Double) Dim w() As Double = New Double(N){} Dim v() As Double = New Double(N){} For k As Integer = 0 To N - 2 d(k) = a(k, k) Dim t As Double = 0.0 For i As Integer = k + 1 To N w(i) = a(i, k) t += w(i) * w(i) Next Dim s As Double = Math.Sqrt(t) If w(k + 1) < 0 Then s = -s End If If Math.Abs(s) < 0.00000000001 Then e(k + 1) = 0.0 Else e(k + 1) = -s w(k + 1) += s s = 1 / Math.Sqrt(w(k + 1) * s) For i As Integer = k + 1 To N w(i) *= s Next For i As Integer = k + 1 To N s = 0.0 For j As Integer = k + 1 To N If j <= i Then s += a(i, j) * w(j) Else s += a(j, i) * w(j) End If Next v(i) = s Next s = 0.0 For i As Integer = k + 1 To N s += w(i) * v(i) Next s /= 2.0 For i As Integer = k + 1 To N v(i) -= s * w(i) Next For i As Integer = k + 1 To N For j As Integer = k + 1 To i a(i, j) -= w(i) * v(j) + w(j) * v(i) Next Next '次の行は固有ベクトルを求めないなら不要 For i As Integer = k + 1 To N a(i, k) = w(i) Next End If Next d(N - 1) = a(N - 1, N - 1) d(N) = a(N, N) e(0) = 0.0 e(N) = a(N, N - 1) '次の行は固有ベクトルを求めないなら不要 For k As Integer = N To 0 Step -1 w = {0.0 ,0.0 ,0.0 ,0.0} If k < N - 1 Then For i As Integer = k + 1 To N w(i) = a(i, k) Next For i As Integer = k + 1 To N Dim s As Double = 0.0 For j As Integer = k + 1 To N s += a(i, j) * w(j) Next v(i) = s Next For i As Integer = k + 1 To N For j As Integer = k + 1 To N a(i, j) -= v(i) * w(j) Next Next End If For i As Integer = 0 To N a(i, k) = 0.0 Next a(k, k) = 1.0 Next End Sub 'QR分解 Private Sub decomp(ByRef a(,) As Double, ByRef d() As Double, ByRef e() As Double) e(0) = 1.0 Dim h As Integer = N Do While (Math.Abs(e(h)) < 0.00000000001) h -= 1 Loop Do While h > 0 e(0) = 0.0 Dim l As Integer = h - 1 Do While (Math.Abs(e(l)) >= 0.00000000001) l -= 1 Loop For j As Integer = 1 To 100 Dim w As Double = (d(h - 1) - d(h)) / 2.0 Dim s As Double = Math.Sqrt(w * w + e(h) * e(h)) If w < 0.0 Then s = -s End If Dim x As Double = d(l) - d(h) + e(h) * e(h) / (w + s) Dim y As Double = e(l + 1) Dim z As Double = 0.0 Dim t As Double = 0.0 Dim u As Double = 0.0 For k As Integer = l To h - 1 If Math.Abs(x) >= Math.Abs(y) Then t = -y / x u = 1 / Math.Sqrt(t * t + 1.0) s = t * u Else t = -x / y s = 1 / Math.Sqrt(t * t + 1.0) If t < 0 Then s = -s End If u = t * s End If w = d(k) - d(k + 1) t = (w * s + 2 * u * e(k + 1)) * s d(k ) = d(k ) - t d(k + 1) = d(k + 1) + t e(k ) = u * e(k) - s * z e(k + 1) = e(k + 1) * (u * u - s * s) + w * s * u '次の行は固有ベクトルを求めないなら不要 For i As Integer = 0 To N x = a(k , i) y = a(k + 1, i) a(k , i) = u * x - s * y a(k + 1, i) = s * x + u * y Next If k < N - 1 Then x = e(k + 1) y = -s * e(k + 2) z = y e(k + 2) = u * e(k + 2) End If Next Console.Write(String.Format("{0,3:D}{1}", j, vbTab)) Call disp_vector(d) '収束判定 If (Math.Abs(e(h)) < 0.00000000001) Then Exit For Next e(0) = 1.0 Do While (Math.Abs(e(h)) < 0.00000000001) h -= 1 Loop Loop '次の行は固有ベクトルを求めないなら不要 For k As Integer = 0 To N - 1 Dim l As Integer = k For i As Integer = k + 1 To N If d(i) > d(l) Then l = i End If Next Dim t As Double = d(k) d(k) = d(l) d(l) = t For i As Integer = 0 To N t = a(k, i) a(k, i) = a(l, i) a(l, i) = t Next Next End Sub '１次元配列を表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '２次元配列を表示 Private Sub disp_matrix(ByVal matrix(,) As Double) For row As Integer = 0 To N For col As Integer = 0 To N Console.Write(String.Format("{0,14:F10}{1}", matrix(row, col), vbTab)) Next Console.WriteLine() Next End Sub End Module
Z:\>vbc -nologo VB1106.vb Z:\>VB1106 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
C#
using System; public class CS1106 { 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[] d = new double[N]; double[] e = new double[N]; // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 decomp(a, d, e); Console.WriteLine(); Console.WriteLine("eigenvalue"); disp_vector(d); Console.WriteLine(); Console.WriteLine("eigenvector"); disp_matrix(a); } // ハウスホルダー変換 private static void tridiagonalize(double[,] a, double[] d, double[] e) { double[] w = new double[N]; double[] v = new double[N]; for (int k = 0; k < N - 2; k++) { d[k] = a[k, k]; double t = 0.0; for (int i = k + 1; i < N; i++) { w[i] = a[i, k]; t += w[i] * w[i]; } double s = Math.Sqrt(t); if (w[k + 1] < 0) s = -s; if (Math.Abs(s) < 0.00000000001) e[k + 1] = 0.0; else { e[k + 1] = -s; w[k + 1] += s; s = 1 / Math.Sqrt(w[k + 1] * s); for (int i = k + 1; i < N; i++) w[i] *= s; for (int i = k + 1; i < N; i++) { s = 0.0; for (int j = k + 1; j < N; j++) { if (j <= i) s += a[i, j] * w[j]; else s += a[j, i] * w[j]; } v[i] = s; } s = 0.0; for (int i = k + 1; i < N; i++) s += w[i] * v[i]; s /= 2.0; for (int i = k + 1; i < N; i++) v[i] -= s * w[i]; for (int i = k + 1; i < N; i++) for (int j = k + 1; j <= i; j++) a[i, j] -= w[i] * v[j] + w[j] * v[i]; // 次の行は固有ベクトルを求めないなら不要 for (int i = k + 1; i < N; i++) a[i, k] = w[i]; } } d[N - 2] = a[N - 2, N - 2]; d[N - 1] = a[N - 1, N - 1]; e[0] = 0.0; e[N - 1] = a[N - 1, N - 2]; // 次の行は固有ベクトルを求めないなら不要 for (int k = N - 1; k >= 0; k--) { w = new double[] {0.0 ,0.0 ,0.0 ,0.0}; if (k < N - 2) { for (int i = k + 1; i < N; i++) w[i] = a[i, k]; for (int i = k + 1; i < N; i++) { double s = 0.0; for (int j = k + 1; j < N; j++) s += a[i, j] * w[j]; v[i] = s; } for (int i = k + 1; i < N; i++) for (int j = k + 1; j < N; j++) a[i, j] -= v[i] * w[j]; } for (int i = 0; i < N; i++) a[i, k] = 0.0; a[k, k] = 1.0; } } // QR分解 private static void decomp(double[,] a, double[] d, double[] e) { e[0] = 1.0; int h = N - 1; while (Math.Abs(e[h]) < 0.00000000001) h--; while (h > 0) { e[0] = 0.0; int l = h - 1; while (Math.Abs(e[l]) >= 0.00000000001) l--; for (int j = 1; j <= 100; j++) { double w = (d[h - 1] - d[h]) / 2.0; double s = Math.Sqrt(w * w + e[h] * e[h]); if (w < 0.0) s = -s; double x = d[l] - d[h] + e[h] * e[h] / (w + s); double y = e[l + 1]; double z = 0.0; double t = 0.0; double u = 0.0; for (int k = l; k < h; k++) { if (Math.Abs(x) >= Math.Abs(y)) { t = -y / x; u = 1 / Math.Sqrt(t * t + 1.0); s = t * u; } else { t = -x / y; s = 1 / Math.Sqrt(t * t + 1.0); if (t < 0) s = -s; u = t * s; } w = d[k] - d[k + 1]; t = (w * s + 2 * u * e[k + 1]) * s; d[k ] = d[k ] - t; d[k + 1] = d[k + 1] + t; e[k ] = u * e[k] - s * z; e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u; // 次の行は固有ベクトルを求めないなら不要 for (int i = 0; i < N; i++) { x = a[k , i]; y = a[k + 1, i]; a[k , i] = u * x - s * y; a[k + 1, i] = s * x + u * y; } if (k < N - 2) { x = e[k + 1]; y = -s * e[k + 2]; z = y; e[k + 2] = u * e[k + 2]; } } Console.Write(string.Format("{0,3:D}\t", j)); disp_vector(d); // 収束判定 if (Math.Abs(e[h]) < 0.00000000001) break; } e[0] = 1.0; while (Math.Abs(e[h]) < 0.00000000001) h--; } // 次の行は固有ベクトルを求めないなら不要 for (int k = 0; k < N - 1; k++) { int l = k; for (int i = k + 1; i < N; i++) if (d[i] > d[l]) l = i; double t = d[k]; d[k] = d[l]; d[l] = t; for (int i = 0; i < N; i++) { t = a[k, i]; a[k, i] = a[l, i]; a[l, i] = t; } } } // １次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // ２次元配列を表示 private static void disp_matrix(double[,] matrix) { for (int row = 0; row < N; row++) { for (int col = 0; col < N; col++) Console.Write(string.Format("{0,14:F10}\t", matrix[row,col])); Console.WriteLine(); } } }
Z:\>csc -nologo CS1106.cs Z:\>CS1106 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Java
import java.lang.*; public class Java1106 { 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[] d = new double[N]; double[] e = new double[N]; // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 decomp(a, d, e); System.out.println(); System.out.println("eigenvalue"); disp_vector(d); System.out.println(); System.out.println("eigenvector"); disp_matrix(a); } // ハウスホルダー変換 private static void tridiagonalize(double[][] a, double[] d, double[] e) { double[] w = new double[N]; double[] v = new double[N]; for (int k = 0; k < N - 2; k++) { d[k] = a[k][k]; double t = 0.0; for (int i = k + 1; i < N; i++) { w[i] = a[i][k]; t += w[i] * w[i]; } double s = Math.sqrt(t); if (w[k + 1] < 0) s = -s; if (Math.abs(s) < 0.00000000001) e[k + 1] = 0.0; else { e[k + 1] = -s; w[k + 1] += s; s = 1 / Math.sqrt(w[k + 1] * s); for (int i = k + 1; i < N; i++) w[i] *= s; for (int i = k + 1; i < N; i++) { s = 0.0; for (int j = k + 1; j < N; j++) { if (j <= i) s += a[i][j] * w[j]; else s += a[j][i] * w[j]; } v[i] = s; } s = 0.0; for (int i = k + 1; i < N; i++) s += w[i] * v[i]; s /= 2.0; for (int i = k + 1; i < N; i++) v[i] -= s * w[i]; for (int i = k + 1; i < N; i++) for (int j = k + 1; j <= i; j++) a[i][j] -= w[i] * v[j] + w[j] * v[i]; // 次の行は固有ベクトルを求めないなら不要 for (int i = k + 1; i < N; i++) a[i][k] = w[i]; } } d[N - 2] = a[N - 2][N - 2]; d[N - 1] = a[N - 1][N - 1]; e[0] = 0.0; e[N - 1] = a[N - 1][N - 2]; // 次の行は固有ベクトルを求めないなら不要 for (int k = N - 1; k >= 0; k--) { w = new double[] {0.0 ,0.0 ,0.0 ,0.0}; if (k < N - 2) { for (int i = k + 1; i < N; i++) w[i] = a[i][k]; for (int i = k + 1; i < N; i++) { double s = 0.0; for (int j = k + 1; j < N; j++) s += a[i][j] * w[j]; v[i] = s; } for (int i = k + 1; i < N; i++) for (int j = k + 1; j < N; j++) a[i][j] -= v[i] * w[j]; } for (int i = 0; i < N; i++) a[i][k] = 0.0; a[k][k] = 1.0; } } // QR分解 private static void decomp(double[][] a, double[] d, double[] e) { e[0] = 1.0; int h = N - 1; while (Math.abs(e[h]) < 0.00000000001) h--; while (h > 0) { e[0] = 0.0; int l = h - 1; while (Math.abs(e[l]) >= 0.00000000001) l--; for (int j = 1; j <= 100; j++) { double w = (d[h - 1] - d[h]) / 2.0; double s = Math.sqrt(w * w + e[h] * e[h]); if (w < 0.0) s = -s; double x = d[l] - d[h] + e[h] * e[h] / (w + s); double y = e[l + 1]; double z = 0.0; double t = 0.0; double u = 0.0; for (int k = l; k < h; k++) { if (Math.abs(x) >= Math.abs(y)) { t = -y / x; u = 1 / Math.sqrt(t * t + 1.0); s = t * u; } else { t = -x / y; s = 1 / Math.sqrt(t * t + 1.0); if (t < 0) s = -s; u = t * s; } w = d[k] - d[k + 1]; t = (w * s + 2 * u * e[k + 1]) * s; d[k ] = d[k ] - t; d[k + 1] = d[k + 1] + t; e[k ] = u * e[k] - s * z; e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u; // 次の行は固有ベクトルを求めないなら不要 for (int i = 0; i < N; i++) { x = a[k ][i]; y = a[k + 1][i]; a[k ][i] = u * x - s * y; a[k + 1][i] = s * x + u * y; } if (k < N - 2) { x = e[k + 1]; y = -s * e[k + 2]; z = y; e[k + 2] = u * e[k + 2]; } } System.out.print(String.format("%3d\t", j)); disp_vector(d); // 収束判定 if (Math.abs(e[h]) < 0.00000000001) break; } e[0] = 1.0; while (Math.abs(e[h]) < 0.00000000001) h--; } // 次の行は固有ベクトルを求めないなら不要 for (int k = 0; k < N - 1; k++) { int l = k; for (int i = k + 1; i < N; i++) if (d[i] > d[l]) l = i; double t = d[k]; d[k] = d[l]; d[l] = t; for (int i = 0; i < N; i++) { t = a[k][i]; a[k][i] = a[l][i]; a[l][i] = t; } } } // １次元配列を表示 private static void disp_vector(double[] row) { for (double col: row) System.out.print(String.format("%14.10f\t", col)); System.out.println(); } // ２次元配列を表示 private static void disp_matrix(double[][] matrix) { for (int row = 0; row < N; row++) { for (int col = 0; col < N; col++) System.out.print(String.format("%14.10f\t", matrix[row][col])); System.out.println(); } } }
Z:\>javac Java1106.java Z:\>java Java1106 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
C++
#include <iostream> #include <iomanip> #include <math.h> using namespace std; const int N = 4; // ハウスホルダー変換 void tridiagonalize(double a[N][N], double d[N], double e[N]); // QR分解 void decomp(double a[N][N], double d[N], double e[N]); // １次元配列を表示 void disp_vector(double row[N]); // ２次元配列を表示 void disp_matrix(double matrix[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 d[N]; double e[N]; // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 decomp(a, d, e); cout << endl << "eigenvalue" << endl; disp_vector(d); cout << endl << "eigenvector" << endl; disp_matrix(a); return 0; } // ハウスホルダー変換 void tridiagonalize(double a[N][N], double d[N], double e[N]) { double v[N]; for (int k = 0; k < N - 2; k++) { double w[N] = {0.0 ,0.0 ,0.0 ,0.0}; d[k] = a[k][k]; double t = 0.0; for (int i = k + 1; i < N; i++) { w[i] = a[i][k]; t += w[i] * w[i]; } double s = sqrt(t); if (w[k + 1] < 0) s = -s; if (fabs(s) < 0.00000000001) e[k + 1] = 0.0; else { e[k + 1] = -s; w[k + 1] += s; s = 1 / sqrt(w[k + 1] * s); for (int i = k + 1; i < N; i++) w[i] *= s; for (int i = k + 1; i < N; i++) { s = 0.0; for (int j = k + 1; j < N; j++) { if (j <= i) s += a[i][j] * w[j]; else s += a[j][i] * w[j]; } v[i] = s; } s = 0.0; for (int i = k + 1; i < N; i++) s += w[i] * v[i]; s /= 2.0; for (int i = k + 1; i < N; i++) v[i] -= s * w[i]; for (int i = k + 1; i < N; i++) for (int j = k + 1; j <= i; j++) a[i][j] -= w[i] * v[j] + w[j] * v[i]; // 次の行は固有ベクトルを求めないなら不要 for (int i = k + 1; i < N; i++) a[i][k] = w[i]; } } d[N - 2] = a[N - 2][N - 2]; d[N - 1] = a[N - 1][N - 1]; e[0] = 0.0; e[N - 1] = a[N - 1][N - 2]; // 次の行は固有ベクトルを求めないなら不要 for (int k = N - 1; k >= 0; k--) { double w[N] = {0.0 ,0.0 ,0.0 ,0.0}; if (k < N - 2) { for (int i = k + 1; i < N; i++) w[i] = a[i][k]; for (int i = k + 1; i < N; i++) { double s = 0.0; for (int j = k + 1; j < N; j++) s += a[i][j] * w[j]; v[i] = s; } for (int i = k + 1; i < N; i++) for (int j = k + 1; j < N; j++) a[i][j] -= v[i] * w[j]; } for (int i = 0; i < N; i++) a[i][k] = 0.0; a[k][k] = 1.0; } } // QR分解 void decomp(double a[N][N], double d[N], double e[N]) { e[0] = 1.0; int h = N - 1; while (fabs(e[h]) < 0.00000000001) h--; while (h > 0) { e[0] = 0.0; int l = h - 1; while (fabs(e[l]) >= 0.00000000001) l--; for (int j = 1; j <= 100; j++) { double w = (d[h - 1] - d[h]) / 2.0; double s = sqrt(w * w + e[h] * e[h]); if (w < 0.0) s = -s; double x = d[l] - d[h] + e[h] * e[h] / (w + s); double y = e[l + 1]; double z = 0.0; double t; double u; for (int k = l; k < h; k++) { if (fabs(x) >= fabs(y)) { t = -y / x; u = 1 / sqrt(t * t + 1.0); s = t * u; } else { t = -x / y; s = 1 / sqrt(t * t + 1.0); if (t < 0) s = -s; u = t * s; } w = d[k] - d[k + 1]; t = (w * s + 2 * u * e[k + 1]) * s; d[k ] = d[k ] - t; d[k + 1] = d[k + 1] + t; e[k ] = u * e[k] - s * z; e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u; // 次の行は固有ベクトルを求めないなら不要 for (int i = 0; i < N; i++) { x = a[k ][i]; y = a[k + 1][i]; a[k ][i] = u * x - s * y; a[k + 1][i] = s * x + u * y; } if (k < N - 2) { x = e[k + 1]; y = -s * e[k + 2]; z = y; e[k + 2] = u * e[k + 2]; } } cout << setw(3) << j << "\t"; disp_vector(d); // 収束判定 if (fabs(e[h]) < 0.00000000001) break; } e[0] = 1.0; while (fabs(e[h]) < 0.00000000001) h--; } // 次の行は固有ベクトルを求めないなら不要 for (int k = 0; k < N - 1; k++) { int l = k; for (int i = k + 1; i < N; i++) if (d[i] > d[l]) l = i; double t = d[k]; d[k] = d[l]; d[l] = t; for (int i = 0; i < N; i++) { t = a[k][i]; a[k][i] = a[l][i]; a[l][i] = t; } } } // １次元配列を表示 void disp_vector(double row[N]) { for (int i = 0; i < N; i++) cout << setw(14) << fixed << setprecision(10) << row[i] << "\t"; cout << endl; } // ２次元配列を表示 void disp_matrix(double matrix[N][N]) { for (int row = 0; row < N; row++) { for (int col = 0; col < N; col++) cout << setw(14) << fixed << setprecision(10) << matrix[row][col] << "\t"; cout << endl; } }
Z:\>bcc32 -q CP1106.cpp cp1106.cpp: Z:\>CP1106 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 -0.0000000000
Objective-C
#import <Foundation/Foundation.h> #import <math.h> #define N 4 // ハウスホルダー変換 void tridiagonalize(double a[N][N], double d[N], double e[N]); // QR分解 void decomp(double a[N][N], double d[N], double e[N]); // １次元配列を表示 void disp_vector(double row[N]); // ２次元配列を表示 void disp_matrix(double matrix[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 d[N]; double e[N]; // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 printf("\nQR分解\n"); decomp(a, d, e); printf("\neigenvalue\n"); disp_vector(d); printf("\neigenvector\n"); disp_matrix(a); return 0; } // ハウスホルダー変換 void tridiagonalize(double a[N][N], double d[N], double e[N]) { int i, j, k; double v[N]; for (k = 0; k < N - 2; k++) { double w[N] = {0.0 ,0.0 ,0.0 ,0.0}; d[k] = a[k][k]; double t = 0.0; for (i = k + 1; i < N; i++) { w[i] = a[i][k]; t += w[i] * w[i]; } double s = sqrt(t); if (w[k + 1] < 0) s = -s; if (fabs(s) < 0.00000000001) e[k + 1] = 0.0; else { e[k + 1] = -s; w[k + 1] += s; s = 1 / sqrt(w[k + 1] * s); for (i = k + 1; i < N; i++) w[i] *= s; for (i = k + 1; i < N; i++) { s = 0.0; for (j = k + 1; j < N; j++) { if (j <= i) s += a[i][j] * w[j]; else s += a[j][i] * w[j]; } v[i] = s; } s = 0.0; for (i = k + 1; i < N; i++) s += w[i] * v[i]; s /= 2.0; for (i = k + 1; i < N; i++) v[i] -= s * w[i]; for (i = k + 1; i < N; i++) for (j = k + 1; j <= i; j++) a[i][j] -= w[i] * v[j] + w[j] * v[i]; // 次の行は固有ベクトルを求めないなら不要 for (i = k + 1; i < N; i++) a[i][k] = w[i]; } } d[N - 2] = a[N - 2][N - 2]; d[N - 1] = a[N - 1][N - 1]; e[0] = 0.0; e[N - 1] = a[N - 1][N - 2]; // 次の行は固有ベクトルを求めないなら不要 for (k = N - 1; k >= 0; k--) { double w[N] = {0.0 ,0.0 ,0.0 ,0.0}; if (k < N - 2) { for (i = k + 1; i < N; i++) w[i] = a[i][k]; for (i = k + 1; i < N; i++) { double s = 0.0; for (j = k + 1; j < N; j++) s += a[i][j] * w[j]; v[i] = s; } for (i = k + 1; i < N; i++) for (j = k + 1; j < N; j++) a[i][j] -= v[i] * w[j]; } for (i = 0; i < N; i++) a[i][k] = 0.0; a[k][k] = 1.0; } } // QR分解 void decomp(double a[N][N], double d[N], double e[N]) { int i, j, k; e[0] = 1.0; int h = N - 1; while (fabs(e[h]) < 0.00000000001) h--; while (h > 0) { e[0] = 0.0; int l = h - 1; while (fabs(e[l]) >= 0.00000000001) l--; for (j = 1; j <= 100; j++) { double w = (d[h - 1] - d[h]) / 2.0; double s = sqrt(w * w + e[h] * e[h]); if (w < 0.0) s = -s; double x = d[l] - d[h] + e[h] * e[h] / (w + s); double y = e[l + 1]; double z = 0.0; double t = 0.0; double u = 0.0; for (k = l; k < h; k++) { if (fabs(x) >= fabs(y)) { t = -y / x; u = 1 / sqrt(t * t + 1.0); s = t * u; } else { t = -x / y; s = 1 / sqrt(t * t + 1.0); if (t < 0) s = -s; u = t * s; } w = d[k] - d[k + 1]; t = (w * s + 2 * u * e[k + 1]) * s; d[k ] = d[k ] - t; d[k + 1] = d[k + 1] + t; e[k ] = u * e[k] - s * z; e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u; // 次の行は固有ベクトルを求めないなら不要 for (i = 0; i < N; i++) { x = a[k ][i]; y = a[k + 1][i]; a[k ][i] = u * x - s * y; a[k + 1][i] = s * x + u * y; } if (k < N - 2) { x = e[k + 1]; y = -s * e[k + 2]; z = y; e[k + 2] = u * e[k + 2]; } } printf("%3d\t", j); disp_vector(d); // 収束判定 if (fabs(e[h]) < 0.00000000001) break; } e[0] = 1.0; while (fabs(e[h]) < 0.00000000001) h--; } // 次の行は固有ベクトルを求めないなら不要 for (k = 0; k < N - 1; k++) { int l = k; for (i = k + 1; i < N; i++) if (d[i] > d[l]) l = i; double t = d[k]; d[k] = d[l]; d[l] = t; for (i = 0; i < N; i++) { t = a[k][i]; a[k][i] = a[l][i]; a[l][i] = t; } } } // １次元配列を表示 void disp_vector(double row[N]) { int i; for (i = 0; i < N; i++) printf("%14.10f\t", row[i]); printf("\n"); } // ２次元配列を表示 void disp_matrix(double matrix[N][N]) { int col, row; for (row = 0; row < N; row++) { for (col = 0; col < N; col++) printf("%14.10f\t", matrix[row][col]); printf("\n"); } }
xxxxxx@yyyyyy /Z $ gcc -o OC1106 OC1106.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC1106 QR分解 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 -0.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 d[N]; double e[N]; // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 decomp(a, d, e); writefln("\neigenvalue"); disp_vector(d); writefln("\neigenvector"); disp_matrix(a); } // ハウスホルダー変換 void tridiagonalize(ref double[N][N] a, ref double[N] d, ref double[N] e) { double v[N]; foreach (k; 0..N - 2) { double w[N] = [0.0 ,0.0 ,0.0 ,0.0]; d[k] = a[k][k]; double t = 0.0; foreach (i; k + 1..N) { w[i] = a[i][k]; t += w[i] * w[i]; } double s = sqrt(t); if (w[k + 1] < 0) s = -s; if (fabs(s) < 0.00000000001) e[k + 1] = 0.0; else { e[k + 1] = -s; w[k + 1] += s; s = 1 / sqrt(w[k + 1] * s); foreach (i; k + 1..N) w[i] *= s; foreach (i; k + 1..N) { s = 0.0; foreach (j; k + 1..N) { if (j <= i) s += a[i][j] * w[j]; else s += a[j][i] * w[j]; } v[i] = s; } s = 0.0; foreach (i; k + 1..N) s += w[i] * v[i]; s /= 2.0; foreach (i; k + 1..N) v[i] -= s * w[i]; foreach (i; k + 1..N) foreach (j; k + 1..i + 1) a[i][j] -= w[i] * v[j] + w[j] * v[i]; // 次の行は固有ベクトルを求めないなら不要 foreach (i; k + 1..N) a[i][k] = w[i]; } } d[N - 2] = a[N - 2][N - 2]; d[N - 1] = a[N - 1][N - 1]; e[0] = 0.0; e[N - 1] = a[N - 1][N - 2]; // 次の行は固有ベクトルを求めないなら不要 foreach_reverse (k; 0..N) { double w[N] = [0.0 ,0.0 ,0.0 ,0.0]; if (k < N - 2) { foreach (i; k + 1..N) w[i] = a[i][k]; foreach (i; k + 1..N) { double s = 0.0; foreach (j; k + 1..N) s += a[i][j] * w[j]; v[i] = s; } foreach (i; k + 1..N) foreach (j; k + 1..N) a[i][j] -= v[i] * w[j]; } foreach (i; 0..N) a[i][k] = 0.0; a[k][k] = 1.0; } } // QR分解 void decomp(ref double[N][N] a, ref double[N] d, ref double[N] e) { e[0] = 1.0; int h = N - 1; while (fabs(e[h]) < 0.00000000001) h--; while (h > 0) { e[0] = 0.0; int l = h - 1; while (fabs(e[l]) >= 0.00000000001) l--; foreach (j; 1..100) { double w = (d[h - 1] - d[h]) / 2.0; double s = sqrt(w * w + e[h] * e[h]); if (w < 0.0) s = -s; double x = d[l] - d[h] + e[h] * e[h] / (w + s); double y = e[l + 1]; double z = 0.0; double t = 0.0; double u = 0.0; foreach (k; l..h) { if (fabs(x) >= fabs(y)) { t = -y / x; u = 1 / sqrt(t * t + 1.0); s = t * u; } else { t = -x / y; s = 1 / sqrt(t * t + 1.0); if (t < 0) s = -s; u = t * s; } w = d[k] - d[k + 1]; t = (w * s + 2 * u * e[k + 1]) * s; d[k ] = d[k ] - t; d[k + 1] = d[k + 1] + t; e[k ] = u * e[k] - s * z; e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u; // 次の行は固有ベクトルを求めないなら不要 foreach (i; 0..N) { x = a[k ][i]; y = a[k + 1][i]; a[k ][i] = u * x - s * y; a[k + 1][i] = s * x + u * y; } if (k < N - 2) { x = e[k + 1]; y = -s * e[k + 2]; z = y; e[k + 2] = u * e[k + 2]; } } writef("%3d\t", j); disp_vector(d); // 収束判定 if (fabs(e[h]) < 0.00000000001) break; } e[0] = 1.0; while (fabs(e[h]) < 0.00000000001) h--; } // 次の行は固有ベクトルを求めないなら不要 foreach (k; 0..N - 1) { int l = k; foreach (i; k + 1..N) if (d[i] > d[l]) l = i; double t = d[k]; d[k] = d[l]; d[l] = t; foreach (i; 0..N) { t = a[k][i]; a[k][i] = a[l][i]; a[l][i] = t; } } } // １次元配列を表示 void disp_vector(double[N] row) { foreach (i; 0..N) writef("%14.10f\t", row[i]); writefln(""); } // ２次元配列を表示 void disp_matrix(double[N][N] matrix) { foreach (row; 0..N) { foreach (col; 0..N) writef("%14.10f\t", matrix[row][col]); writefln(""); } }
Z:\>dmd D1106.d Z:\>D1106 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 -0.0000000000 -0.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 d = []float64 {0.0 ,0.0 ,0.0 ,0.0} var e = []float64 {0.0 ,0.0 ,0.0 ,0.0} // ハウスホルダー変換 tridiagonalize(&a, d, e) // QR分解 decomp(&a, d, e) fmt.Println("\neigenvalue") disp_vector(d) fmt.Println("\neigenvector") disp_matrix(a) } // ハウスホルダー変換 func tridiagonalize(a *[N][N]float64, d []float64, e []float64) { var v = []float64 {0.0 ,0.0 ,0.0 ,0.0} for k := 0; k < N - 2; k++ { var w = []float64 {0.0 ,0.0 ,0.0 ,0.0} d[k] = a[k][k] var t float64 = 0.0 for i := k + 1; i < N; i++ { w[i] = a[i][k] t += w[i] * w[i] } var s float64 = math.Sqrt(t) if (w[k + 1] < 0) { s = -s } if (math.Abs(s) < 0.00000000001) { e[k + 1] = 0.0 } else { e[k + 1] = -s w[k + 1] += s s = 1 / math.Sqrt(w[k + 1] * s) for i := k + 1; i < N; i++ { w[i] *= s } for i := k + 1; i < N; i++ { s = 0.0 for j := k + 1; j < N; j++ { if (j <= i) { s += a[i][j] * w[j] } else { s += a[j][i] * w[j] } } v[i] = s } s = 0.0 for i := k + 1; i < N; i++ { s += w[i] * v[i] } s /= 2.0 for i := k + 1; i < N; i++ { v[i] -= s * w[i] } for i := k + 1; i < N; i++ { for j := k + 1; j < i + 1; j++ { a[i][j] -= w[i] * v[j] + w[j] * v[i] } } // 次の行は固有ベクトルを求めないなら不要 for i := k + 1; i < N; i++ { a[i][k] = w[i] } } } d[N - 2] = a[N - 2][N - 2] d[N - 1] = a[N - 1][N - 1] e[0] = 0.0 e[N - 1] = a[N - 1][N - 2] // 次の行は固有ベクトルを求めないなら不要 for k := N - 1; k >= 0; k-- { var w = []float64 {0.0 ,0.0 ,0.0 ,0.0} if (k < N - 2) { for i := k + 1; i < N; i++ { w[i] = a[i][k] } for i := k + 1; i < N; i++ { var s float64 = 0.0 for j := k + 1; j < N; j++ { s += a[i][j] * w[j] } v[i] = s } for i := k + 1; i < N; i++ { for j := k + 1; j < N; j++ { a[i][j] -= v[i] * w[j] } } } for i := 0; i < N; i++ { a[i][k] = 0.0 } a[k][k] = 1.0 } } // QR分解 func decomp(a *[N][N]float64, d []float64, e []float64) { e[0] = 1.0 var h int = N - 1 for (math.Abs(e[h]) < 0.00000000001) { h-- } for (h > 0) { e[0] = 0.0 var l int = h - 1 for (math.Abs(e[l]) >= 0.00000000001) { l-- } for j := 1; j < 100; j++ { var w float64 = (d[h - 1] - d[h]) / 2.0 var s float64 = math.Sqrt(w * w + e[h] * e[h]) if (w < 0.0) { s = -s } var x float64 = d[l] - d[h] + e[h] * e[h] / (w + s) var y float64 = e[l + 1] var z float64 = 0.0 var t float64 = 0.0 var u float64 = 0.0 for k := l; k < h; k++ { if (math.Abs(x) >= math.Abs(y)) { t = -y / x u = 1 / math.Sqrt(t * t + 1.0) s = t * u } else { t = -x / y s = 1 / math.Sqrt(t * t + 1.0) if (t < 0) { s = -s } u = t * s } w = d[k] - d[k + 1] t = (w * s + 2 * u * e[k + 1]) * s d[k ] = d[k ] - t d[k + 1] = d[k + 1] + t e[k ] = u * e[k] - s * z e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for i := 0; i < N; i++ { x = a[k ][i] y = a[k + 1][i] a[k ][i] = u * x - s * y a[k + 1][i] = s * x + u * y } if (k < N - 2) { x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] } } fmt.Printf("%3d\t", j) disp_vector(d) // 収束判定 if (math.Abs(e[h]) < 0.00000000001) { break } } e[0] = 1.0 for (math.Abs(e[h]) < 0.00000000001) { h-- } } // 次の行は固有ベクトルを求めないなら不要 for k := 0; k < N - 1; k++ { var l int = k for i := k + 1; i < N; i++ { if (d[i] > d[l]) { l = i } } var t float64 = d[k] d[k] = d[l] d[l] = t for i := 0; i < N; i++ { t = a[k][i] a[k][i] = a[l][i] a[l][i] = t } } } // １次元配列を表示 func disp_vector(row[]float64) { for _, col := range row { fmt.Printf("%14.10f\t", col) } fmt.Println("") } // ２次元配列を表示 func disp_matrix(matrix[N][N]float64) { for row := 0; row < N; row++ { for col := 0; col < N; col++ { fmt.Printf("%14.10f\t", matrix[row][col]) } fmt.Println() } }
Z:\>go run GO1106.go 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
Scala
object Scala1106 { 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 (d, e) = tridiagonalize(a) // QR分解 decomp(a, d, e) println() println("eigenvalue") disp_vector(d) println() println("eigenvector") disp_matrix(a) } // ハウスホルダー変換 def tridiagonalize(a: => Array[Array[Double]]): (Array[Double], Array[Double]) = { var v = new Array[Double](N + 1) var d = new Array[Double](N + 1) var e = new Array[Double](N + 1) for (k <- 0 to N - 2) { d(k) = a(k)(k) var w = (for (i <- 0 to N) yield (if (i < k + 1) 0.0 else a(i)(k)) ).toArray var s = Math.sqrt(w.map(n => n * n).sum) if (w(k + 1) < 0) s = -s if (Math.abs(s) < 0.00000000001) e(k + 1) = 0.0 else { e(k + 1) = -s w(k + 1) += s s = 1 / Math.sqrt(w(k + 1) * s) for (i <- k + 1 to N) w(i) *= s for (i <- k + 1 to N) { v(i) = (for (j <- k + 1 to N) yield (if (j <= i) a(i)(j) * w(j) else a(j)(i) * w(j)) ).sum } s = (for (i <- k + 1 to N) yield (w(i) * v(i)) ).sum / 2.0 for (i <- k + 1 to N) v(i) -= s * w(i) for (i <- k + 1 to N) for (j <- k + 1 to i) a(i)(j) -= w(i) * v(j) + w(j) * v(i) // 次の行は固有ベクトルを求めないなら不要 for (i <- k + 1 to N) a(i)(k) = w(i) } } d(N - 1) = a(N - 1)(N - 1) d(N ) = a(N )(N ) e(0) = 0.0 e(N ) = a(N )(N - 1) // 次の行は固有ベクトルを求めないなら不要 for (k <- (N to 0 by -1)) { var w = (for (i <- 0 to N) yield (if (i < k + 1) 0.0 else a(i)(k)) ).toArray for (i <- k + 1 to N) { v(i) = (for (j <- k + 1 to N) yield (a(i)(j) * w(j)) ).sum } for (i <- k + 1 to N) for (j <- k + 1 to N) a(i)(j) -= v(i) * w(j) for (i <- 0 to N) a(i)(k) = 0.0 a(k)(k) = 1.0 } return (d, e) } // QR分解 def decomp(a: => Array[Array[Double]], d: => Array[Double], e: => Array[Double]) { e(0) = 1.0 var h = N while (Math.abs(e(h)) < 0.00000000001) h -= 1 while (h > 0) { e(0) = 0.0 var l = h - 1 while (Math.abs(e(l)) >= 0.00000000001) l -= 1 var j = 1 do { var w = (d(h - 1) - d(h)) / 2.0 var s = Math.sqrt(w * w + e(h) * e(h)) if (w < 0.0) s = -s var x = d(l) - d(h) + e(h) * e(h) / (w + s) var y = e(l + 1) var z = 0.0 var t = 0.0 var u = 0.0 for (k <- l to h - 1) { if (Math.abs(x) >= Math.abs(y)) { t = -y / x u = 1 / Math.sqrt(t * t + 1.0) s = t * u } else { t = -x / y s = 1 / Math.sqrt(t * t + 1.0) * (if (t < 0) -1.0 else 1.0) u = t * s } w = d(k) - d(k + 1) t = (w * s + 2 * u * e(k + 1)) * s d(k ) = d(k ) - t d(k + 1) = d(k + 1) + t e(k ) = u * e(k) - s * z e(k + 1) = e(k + 1) * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for (i <- 0 to N) { x = a(k )(i) y = a(k + 1)(i) a(k )(i) = u * x - s * y a(k + 1)(i) = s * x + u * y } if (k < N - 2) { x = e(k + 1) y = -s * e(k + 2) z = y e(k + 2) = u * e(k + 2) } } print("%3d\t".format(j)) j += 1 disp_vector(d) } while (j <= 100 && Math.abs(e(h)) >= 0.00000000001) e(0) = 1.0 while (Math.abs(e(h)) < 0.00000000001) h -= 1 } // 次の行は固有ベクトルを求めないなら不要 for (k <- 0 to N - 1) { var l = k for (i <- k + 1 to N) if (d(i) > d(l)) l = i var t = d(k) d(k) = d(l) d(l) = t for (i <- 0 to N) { t = a(k)(i) a(k)(i) = a(l)(i) a(l)(i) = t } } } // １次元配列を表示 def disp_vector(row:Array[Double]) = { row.map( col => print("%14.10f\t".format(col))) println() } // ２次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) { for (row <- 0 to N) { for (col <- 0 to N) print("%14.10f\t".format(matrix(row)(col))) println() } } }
Z:\>scala Scala1106.scala 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 -0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000
F#
module Fs1106 open System let N = 3 // １次元配列を表示 let disp_vector (vector:float[]) = vector |> Array.iter (fun col -> printf "%14.10f" col) printfn "" // ２次元配列を表示 let disp_matrix (matrix:float[][]) = for row in [0..N] do for col in [0..N] do printf "%14.10f" matrix.[row].[col] printfn "" // ハウスホルダー変換 let tridiagonalize (a:float[][]) : float[] * float[] = let v:float[] = [| for i in 0..N -> 0.0 |] let d:float[] = [| for i in 0..N -> 0.0 |] let e:float[] = [| for i in 0..N -> 0.0 |] for k in [0..(N - 2)] do d.[k] <- a.[k].[k] let w = [| for i in [0..N] -> if (i < k + 1) then 0.0 else a.[i].[k] |] let mutable s = Math.Sqrt(w |> Array.map (fun n -> n * n) |> Array.sum) if (w.[k + 1] < 0.0) then s <- -s if (Math.Abs(s) < 0.00000000001) then e.[k + 1] <- 0.0 else e.[k + 1] <- -s w.[k + 1] <- w.[k + 1] + s s <- (1.0 / Math.Sqrt(w.[k + 1] * s)) for i in [(k + 1)..N] do w.[i] <- w.[i] * s for i in [(k + 1)..N] do v.[i] <- [| for j in [(k + 1)..N] -> if (j <= i) then a.[i].[j] * w.[j] else a.[j].[i] * w.[j] |] |> Array.sum s <- ( [| for i in[(k + 1)..N] -> w.[i] * v.[i] |] |> Array.sum ) / 2.0 for i in [(k + 1)..N] do v.[i] <- v.[i] - s * w.[i] for i in [(k + 1)..N] do for j in [(k + 1)..i] do a.[i].[j] <- a.[i].[j] - (w.[i] * v.[j] + w.[j] * v.[i]) // 次の行は固有ベクトルを求めないなら不要 for i in [(k + 1)..N] do a.[i].[k] <- w.[i] d.[N - 1] <- a.[N - 1].[N - 1] d.[N ] <- a.[N ].[N ] e.[0] <- 0.0 e.[N ] <- a.[N ].[N - 1] // 次の行は固有ベクトルを求めないなら不要 for k = N downto 0 do if k < N - 1 then let w = [| for i in [0..N] -> if (i < k + 1) then 0.0 else a.[i].[k] |] for i in [(k + 1)..N] do v.[i] <- [| for j in [(k + 1)..N] -> a.[i].[j] * w.[j] |] |> Array.sum for i in [(k + 1)..N] do for j in [(k + 1)..N] do a.[i].[j] <- a.[i].[j] - v.[i] * w.[j] for i in [0..N] do a.[i].[k] <- 0.0 a.[k].[k] <- 1.0 (d, e) // QR分解 let decomp (a:float[][]) (d:float[]) (e:float[]) = e.[0] <- 1.0 let mutable h = N while (Math.Abs(e.[h]) < 0.00000000001) do h <- h - 1 while (h > 0) do e.[0] <- 0.0 let mutable l = h - 1 while (Math.Abs(e.[l]) >= 0.00000000001) do l <- l - 1 let mutable j = 1 let mutable sw = true while sw || (j <= 100 && Math.Abs(e.[h]) >= 0.00000000001) do if sw then sw <- false let mutable w = (d.[h - 1] - d.[h]) / 2.0 let mutable s = Math.Sqrt(w * w + e.[h] * e.[h]) if (w < 0.0) then s <- -s let mutable x = d.[l] - d.[h] + e.[h] * e.[h] / (w + s) let mutable y = e.[l + 1] let mutable z = 0.0 let mutable t = 0.0 let mutable u = 0.0 for k in [l..(h - 1)] do if (Math.Abs(x) >= Math.Abs(y)) then t <- -y / x u <- 1.0 / Math.Sqrt(t * t + 1.0) s <- t * u else t <- -x / y s <- 1.0 / Math.Sqrt(t * t + 1.0) * if (t < 0.0) then -1.0 else 1.0 u <- t * s w <- d.[k] - d.[k + 1] t <- (w * s + 2.0 * u * e.[k + 1]) * s d.[k ] <- d.[k ] - t d.[k + 1] <- d.[k + 1] + t e.[k ] <- u * e.[k] - s * z e.[k + 1] <- e.[k + 1] * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for i in [0..N] do x <- a.[k ].[i] y <- a.[k + 1].[i] a.[k ].[i] <- u * x - s * y a.[k + 1].[i] <- s * x + u * y if k < N - 2 then x <- e.[k + 1] y <- -s * e.[k + 2] z <- y e.[k + 2] <- u * e.[k + 2] printf "%3d\t" j j <- j + 1 disp_vector d e.[0] <- 1.0 while (Math.Abs(e.[h]) < 0.00000000001) do h <- h - 1 // 次の行は固有ベクトルを求めないなら不要 for k in [0..(N - 1)] do let mutable l = k for i in [(k + 1)..N] do if d.[i] > d.[l] then l <- i let t = d.[k] d.[k] <- d.[l] d.[l] <- t for i in [0..N] do let t = a.[k].[i] a.[k].[i] <- a.[l].[i] a.[l].[i] <- t // ハウスホルダー変換と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 (d, e) = tridiagonalize a // QR分解 decomp a d e printfn "" printfn "eigenvalue" disp_vector d printfn "" printfn "eigenvector" disp_matrix a exit 0
Z:\>fsi --nologo --quiet Fs1106.fs 1 7.8421052632 3.1783029001 4.9795918367 2.0000000000 2 8.4366343427 2.5633699270 4.9999957303 2.0000000000 3 8.9326997461 2.0673002539 5.0000000000 2.0000000000 4 9.2864656206 1.7135343794 5.0000000000 2.0000000000 1 10.0000000000 1.0000000000 5.0000000000 2.0000000000 eigenvalue 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 0.6324555320 0.6324555320 0.3162277660 0.3162277660 0.3162277660 0.3162277660 -0.6324555320 -0.6324555320 0.0000000000 0.0000000000 -0.7071067812 0.7071067812 0.7071067812 -0.7071067812 0.0000000000 0.0000000000