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;
}
}
}
# 2次元配列を表示
sub disp_matrix
{
my ($matrix) = @_;
foreach $row (@$matrix)
{
foreach $col (@$row)
{
printf("%14.10f\t", $col);
}
print "\n";
}
}
# 1次元配列を表示
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;
}
}
}
# 2次元配列を表示
sub disp_matrix
{
my ($matrix) = @_;
foreach $row (@$matrix)
{
foreach $col (@$row)
{
printf("%14.10f\t", $col);
}
print "\n";
}
}
# 1次元配列を表示
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);
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;
}
}
}
# 2次元配列を表示
sub disp_matrix
{
my ($matrix) = @_;
foreach $row (@$matrix)
{
foreach $col (@$row)
{
printf("%14.10f\t", $col);
}
print "\n";
}
}
# 1次元配列を表示
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];
}
}
# 2次元配列を表示
sub disp_matrix
{
my ($matrix) = @_;
foreach $row (@$matrix)
{
foreach $col (@$row)
{
printf("%14.10f\t", $col);
}
print "\n";
}
}
# 1次元配列を表示
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];
}
}
# 2次元配列を表示
sub disp_matrix
{
my ($matrix) = @_;
foreach $row (@$matrix)
{
foreach $col (@$row)
{
printf("%14.10f\t", $col);
}
print "\n";
}
}
# 1次元配列を表示
sub disp_vector
{
my ($row) = @_;
foreach $col (@$row)
{
printf("%14.10f\t", $col);
}
print "\n";
}