Skip to content

Instantly share code, notes, and snippets.

@fornext1119

fornext1119/index.html

Last active January 14, 2016 00:44
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save fornext1119/14f6d0f26259b7b50643 to your computer and use it in GitHub Desktop.
Save fornext1119/14f6d0f26259b7b50643 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile1.md

ヤコビの反復法

use constant N => 3;

my @a = ([9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6]);
my @b = (20,16,8,17);
my @c = (0,0,0,0);

# ヤコビ反復法
jacobi(\@a, \@b, \@c);

print "解\n";
for $i (0..N)
{
    printf("%14.10f\t", $c[$i]);
}
print "\n";

# ヤコビの反復法
sub jacobi
{
    my ($a, $b, $x0) = @_;

    while (1)
    {
        my @x1 = ();
        my $finish = 1;
        for $i (0..N)
        {
            $x1[$i] = 0;
            for $j (0..N)
            {
                if ($j != $i)
                {
                    $x1[$i] += $$a[$i][$j] * $$x0[$j];
                }
            }
            $x1[$i] = ($$b[$i] - $x1[$i]) / $$a[$i][$i];
            $finish = 0 if (abs($x1[$i] - $$x0[$i]) > 0.0000000001)
        }

        for $i (0..N)
        {
            $$x0[$i] = $x1[$i];
            printf("%14.10f\t", $$x0[$i]);
        }
        print "\n";
        return if ($finish);
    }
}

ガウス・ザイデル法

use constant N => 3;

my @a = ([9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6]);
my @b = (20,16,8,17);
my @c = (0,0,0,0);

# ガウス・ザイデル法
gauss(\@a, \@b, \@c);

print "解\n";
for $i (0..N)
{
    printf("%14.10f\t", $c[$i]);
}
print "\n";

# ガウス・ザイデル法
sub gauss
{
    my ($a, $b, $x0) = @_;

    while (1)
    {
        my $finish = 1;
        for $i (0..N)
        {
            my $x1 = 0;
            for $j (0..N)
            {
                if ($j != $i)
                {
                    $x1 += $$a[$i][$j] * $$x0[$j];
                }
            }
            $x1 = ($$b[$i] - $x1) / $$a[$i][$i];
            $finish = 0 if (abs($x1 - $$x0[$i]) > 0.0000000001);
            $$x0[$i] = $x1;
        }

        for $i (0..N)
        {
            printf("%14.10f\t", $$x0[$i]);
        }
        print "\n";
        return if ($finish);
    }
}

ガウスの消去法

use constant N => 3;

my @a = ([-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]);
my @b = (8,17,20,16);

# ピボット選択
pivoting(\@a,\@b);
print "ピボット選択後\n";
print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# 前進消去
forward_elimination(\@a,\@b);
print "前進消去後\n";
print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# 後退代入
backward_substitution(\@a,\@b);
print "X\n";
disp_vector(\@b);

# 前進消去
sub forward_elimination
{
    my ($a, $b) = @_;

    for $pivot (0..(N - 1))
    {
        for $row (($pivot + 1)..N)
        {
            my $s = $$a[$row][$pivot] / $$a[$pivot][$pivot];
            for $col ($pivot..N)
            {
                $$a[$row][$col] -= $$a[$pivot][$col]	* $s;
            }
            $$b[$row]           -= $$b[$pivot]         	* $s;
        }
    }
}
# 後退代入
sub backward_substitution
{
    my ($a, $b) = @_;

    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $$b[$row] -= $$a[$row][$col] * $$b[$col];
        }
        $$b[$row] /= $$a[$row][$row];
    }
}
# ピボット選択
sub pivoting
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        # 各列で 一番値が大きい行を 探す
        my $max_row =   $pivot;
        my $max_val =   0;
        for $row ($pivot..N)
        {
            if (abs($$a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($$a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            my $tmp;
            for $col (0..N)
            {
                $tmp                =   $$a[$max_row][$col];
                $$a[$max_row][$col] =   $$a[$pivot][$col];
                $$a[$pivot][$col]   =   $tmp;
            }
            $tmp          =   $$b[$max_row];
            $$b[$max_row] =   $$b[$pivot];
            $$b[$pivot]   =   $tmp;
        }
    }
}
# ２次元配列を表示
sub disp_matrix
{
	my ($matrix) = @_;
	foreach $row (@$matrix)
	{ 
		foreach $col (@$row)
		{
			printf("%14.10f\t", $col);
		}    
		print "\n";
	}
}
# １次元配列を表示
sub disp_vector
{
	my ($row) = @_;
	foreach $col (@$row)
	{
		printf("%14.10f\t", $col);
	}    
	print "\n";
}

ガウス・ジョルダン法

use constant N => 3;

my @a = ([-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]);
my @b = (8,17,20,16);

# ピボット選択
pivoting(\@a,\@b);
print "ピボット選択後\n";
print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# 前進消去
forward_elimination(\@a,\@b);
print "前進消去後\n";
print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# 後退代入
backward_substitution(\@a,\@b);
print "X\n";
disp_vector(\@b);

# 前進消去
sub forward_elimination
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        for $row (0..N)
        {
            next if ($row == $pivot);

            my $s = $$a[$row][$pivot] / $$a[$pivot][$pivot];
            for $col ($pivot..N)
            {
                $$a[$row][$col] -= $$a[$pivot][$col] * $s;
            }
            $$b[$row] -= $$b[$pivot] * $s;
        }
    }
}
# 後退代入
sub backward_substitution
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        $$b[$pivot]  /= $$a[$pivot][$pivot];
    }
}

# ピボット選択
sub pivoting
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        # 各列で 一番値が大きい行を 探す
        my $max_row =   $pivot;
        my $max_val =   0;
        for $row ($pivot..N)
        {
            if (abs($$a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($$a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            my $tmp;
            for $col (0..N)
            {
                $tmp                =   $$a[$max_row][$col];
                $$a[$max_row][$col] =   $$a[$pivot][$col];
                $$a[$pivot][$col]   =   $tmp;
            }
            $tmp          =   $$b[$max_row];
            $$b[$max_row] =   $$b[$pivot];
            $$b[$pivot]   =   $tmp;
        }
    }
}
# ２次元配列を表示
sub disp_matrix
{
	my ($matrix) = @_;
	foreach $row (@$matrix)
	{ 
		foreach $col (@$row)
		{
			printf("%14.10f\t", $col);
		}    
		print "\n";
	}
}
# １次元配列を表示
sub disp_vector
{
	my ($row) = @_;
	foreach $col (@$row)
	{
		printf("%14.10f\t", $col);
	}    
	print "\n";
}

LU分解

use constant N => 3;

my @a = ([-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]); 
my @b = (8,17,20,16);

my @a = ([-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]);
my @b = (8,17,20,16);

# ピボット選択
pivoting(\@a,\@b);
print "ピボット選択後\n";
print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# LU分解
# 前進消去
forward_elimination(\@a,\@b);
print "LU\n";
disp_matrix(\@a);

# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);
print "Y\n";
disp_vector(\@y);

# Ux=y から x を求める (後退代入)
my @x = (0,0,0,0);
backward_substitution(\@a,\@y,\@x);
print "X\n";
disp_vector(\@x);

# 前進消去
sub forward_elimination
{
    my ($a, $b) = @_;

    for $pivot (0..(N - 1))
    {
        for $row (($pivot + 1)..N)
        {
            my $s = $$a[$row][$pivot] / $$a[$pivot][$pivot];
            for $col ($pivot..N)
            {
                $$a[$row][$col] -= $$a[$pivot][$col] * $s; # これが 上三角行列
            }
            $$a[$row][$pivot]  = $s;                       # これが 下三角行列
            # $$b[$row]       -= $$b[$pivot] * $s;         # この値は変更しない
        }
    }
}
# Ly=b から y を求める (前進代入)
sub forward_substitution
{
    my ($a, $b, $y) = @_;

    for $row (0..N)
    {
        for $col (0..$row)
        {
            $$b[$row] -= $$a[$row][$col] * $$y[$col];
        }
        $$y[$row] = $$b[$row];
    }
}
# Ux=y から x を求める (後退代入)
sub backward_substitution
{
    my ($a, $y, $x) = @_;

    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $$y[$row] -= $$a[$row][$col] * $$x[$col];
        }
        $$x[$row] = $$y[$row] / $$a[$row][$row];
    }
}

# ピボット選択
sub pivoting
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        # 各列で 一番値が大きい行を 探す
        my $max_row =   $pivot;
        my $max_val =   0;
        for $row ($pivot..N)
        {
            if (abs($$a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($$a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            my $tmp;
            for $col (0..N)
            {
                $tmp                =   $$a[$max_row][$col];
                $$a[$max_row][$col] =   $$a[$pivot][$col];
                $$a[$pivot][$col]   =   $tmp;
            }
            $tmp          =   $$b[$max_row];
            $$b[$max_row] =   $$b[$pivot];
            $$b[$pivot]   =   $tmp;
        }
    }
}
# ２次元配列を表示
sub disp_matrix
{
	my ($matrix) = @_;
	foreach $row (@$matrix)
	{ 
		foreach $col (@$row)
		{
			printf("%14.10f\t", $col);
		}    
		print "\n";
	}
}
# １次元配列を表示
sub disp_vector
{
	my ($row) = @_;
	foreach $col (@$row)
	{
		printf("%14.10f\t", $col);
	}    
	print "\n";
}

コレスキー法

use constant N => 3;

my @a = ([5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]); 
my @b = (34,68,96,125);

print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# LL^T分解
# 前進消去
forward_elimination(\@a,\@b);
print "LL^T\n";
disp_matrix(\@a);

# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);
print "Y\n";
disp_vector(\@y);

# Ux=y から x を求める (後退代入)
my @x = (0,0,0,0);
backward_substitution(\@a,\@y,\@x);
print "X\n";
disp_vector(\@x);

# 前進消去
sub forward_elimination
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        my $s = 0;
        for $col (0..($pivot-1))
        {
            $s += $$a[$pivot][$col] * $$a[$pivot][$col];
        }
        # ここで根号の中が負の値になると計算できない！
        $$a[$pivot][$pivot] = sqrt($$a[$pivot][$pivot] - $s);

        for $row (($pivot + 1)..N)
        {
            $s = 0;
            for $col (0..($pivot-1))
            {
                $s += $$a[$row][$col] * $$a[$pivot][$col];
            }
            $$a[$row][$pivot] = ($$a[$row][$pivot] - $s) / $$a[$pivot][$pivot];
            $$a[$pivot][$row] = $$a[$row][$pivot];
        } 
    }
}
# Ly=b から y を求める (前進代入)
sub forward_substitution
{
    my ($a, $b, $y) = @_;

    for $row (0..N)
    {
        for $col (0..$row)
        {
            $$b[$row] -= $$a[$row][$col] * $$y[$col];
        }
        $$y[$row] = $$b[$row] / $$a[$row][$row];
    }
}
# Ux=y から x を求める (後退代入)
sub backward_substitution
{
    my ($a, $y, $x) = @_;

    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $$y[$row] -= $$a[$row][$col] * $$x[$col];
        }
        $$x[$row] = $$y[$row] / $$a[$row][$row];
    }
}
# ２次元配列を表示
sub disp_matrix
{
	my ($matrix) = @_;
	foreach $row (@$matrix)
	{ 
		foreach $col (@$row)
		{
			printf("%14.10f\t", $col);
		}    
		print "\n";
	}
}
# １次元配列を表示
sub disp_vector
{
	my ($row) = @_;
	foreach $col (@$row)
	{
		printf("%14.10f\t", $col);
	}    
	print "\n";
}

修正コレスキー法

use constant N => 3;

my @a = ([5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]); 
my @b = (34,68,96,125);

print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# LDL^T分解
# 前進消去
forward_elimination(\@a,\@b);
print "LDL^T\n";
disp_matrix(\@a);

# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);
print "Y\n";
disp_vector(\@y);

# Ux=y から x を求める (後退代入)
my @x = (0,0,0,0);
backward_substitution(\@a,\@y,\@x);
print "X\n";
disp_vector(\@x);

# 前進消去
sub forward_elimination
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        # pivot < k の場合
        my $s = 0;
        for $col (0..($pivot-1))
        {
            $s = $$a[$pivot][$col];
            for $k (0..($col-1))
            {
                $s -= $$a[$pivot][$k] * $$a[$col][$k] * $$a[$k][$k];
            }
            $$a[$pivot][$col] = $s / $$a[$col][$col];
            $$a[$col][$pivot] = $$a[$pivot][$col];
        }
 
        # pivot == k の場合
        $s = $$a[$pivot][$pivot];
        for $k (0..($pivot-1))
        {
            $s -= $$a[$pivot][$k] * $$a[$pivot][$k] * $$a[$k][$k];
        }
        $$a[$pivot][$pivot] = $s;
    }
}
# Ly=b から y を求める (前進代入)
sub forward_substitution
{
    my ($a, $b, $y) = @_;

    for $row (0..N)
    {
        for $col (0..$row)
        {
            $$b[$row] -= $$a[$row][$col] * $$y[$col];
        }
        $$y[$row] = $$b[$row];
    }
}
# Ux=y から x を求める (後退代入)
sub backward_substitution
{
    my ($a, $y, $x) = @_;

    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $$y[$row] -= $$a[$row][$col] * $$a[$row][$row] * $$x[$col];
        }
        $$x[$row] = $$y[$row] / $$a[$row][$row];
    }
}
# ２次元配列を表示
sub disp_matrix
{
	my ($matrix) = @_;
	foreach $row (@$matrix)
	{ 
		foreach $col (@$row)
		{
			printf("%14.10f\t", $col);
		}    
		print "\n";
	}
}
# １次元配列を表示
sub disp_vector
{
	my ($row) = @_;
	foreach $col (@$row)
	{
		printf("%14.10f\t", $col);
	}    
	print "\n";
}
Raw
index.html
<html>
<body>
<h1>
hello world!
</h1>
</body>
</html>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment