home > さまざまな言語で数値計算 > 固有値・固有ベクトル >

さまざまな言語で数値計算

Only Do What Only You Can Do

ヤコビ法

ヤコビ法を使って, 全ての固有値と, 対応する固有ベクトルを求める

対称行列 $ \boldsymbol{A} $ の非対角成分のうち絶対値が最大の成分を $ a_{ij} $ とするとき,

その他の対角成分を $ 1 $ , その他の非対角成分を $ 0 $ とした行列 $ \boldsymbol{P} $ を考える.

$ \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

Clojure

Haskell

inserted by FC2 system