print 3 + 5, "\n";
print 3 - 5, "\n";
print 3 * 5, "\n";
print 3 ** 5, "\n";
print 5 / 3, "\n";
print int(5 / 3), "\n";
print 5 % 3, "\n";
printf "%3d\n", 3 * 5;
printf "%23.20f\n", 5 / 3;
my $i = 3 * 5;
print "3 * 5 = ".$i."\n";
print "3 * 5 = ", $i, "\n";
print "3 * 5 = $i\n";
printf "3 * 5 = %d\n", $i;
for (1..9)
{
print "$_, ";
}
print "\n";
for (1..9)
{
if ($_ % 3 == 0)
{
print "$_, ";
}
}
print "\n";
for (1..9)
{
print "$_, " if ($_ % 3 == 0)
}
print "\n";
for (1..9)
{
print "$_, " unless ($_ % 3)
}
print "\n";
use strict;
use warnings;
my $sum = 0;
for (1..99)
{
if ($_ % 3 == 0)
{
$sum += $_;
}
}
print "$sum\n";
map {print $_, ","} 1..9
map {print $_, ","} grep {$_ % 3 == 0} 1..9
use List::Util;
print List::Util::reduce {$a + $b} grep {$_ % 3 == 0} 1..99
print List::Util::sum grep {$_ % 3 == 0} 1..99
use strict;
use warnings;
# 初項:a, 公差:a で 上限:lim の数列の総和を返す関数
sub sn
{
my ($a, $lim) = @_;
my $n = int($lim / $a); # 項数:n = 上限:lim / 公差:a
my $l = $n * $a; # 末項:l = 項数:n * 公差:a
($a + $l) * $n / 2; # 総和:sn = (初項:a + 末項:l) * 項数:n / 2
}
# 3 の倍数の合計
my $sum = sn(3, 999);
print $sum, "\n";
# 10000 までの 自然数の和
# 項目数 n = 10000
my $n = 10000;
print($n * ($n + 1) / 2, "\n");
# 10000 までの 偶数の和
# 項目数 n = 5000
my $n = 10000 / 2;
print($n * ($n + 1), "\n");
# 10000 までの 奇数の和
# 項目数 n = 5000
my $n = 10000 / 2;
print($n ** 2, "\n");
# 1000 までの 自然数の２乗の和
my $n = 1000;
print($n * ($n + 1) * (2 * $n + 1) / 6, "\n");
# 100 までの 自然数の３乗の和
my $n = 100;
print(($n ** 2) * (($n + 1) ** 2) / 4, "\n");
# 初項 2, 公比 3, 項数 10 の等比数列の和
my $n = 10;
my $a = 2;
my $r = 3;
print(($a * (($r ** $n) - 1)) / ($r - 1), "\n");
$a = 5; # 初項 5
$d = 3; # 公差 3
$n = 10; # 項数 10
$p = 1; # 積
for (1..$n)
{
$m = $a + ($d * ($_ - 1));
$p *= $m;
}
print $p;
sub product
{
my ($m, $d, $n) = @_;
if ($n == 0)
{
1
}
else
{
$m * product($m + $d, $d, $n - 1)
}
}
# 初項 5, 公差 3, 項数 10 の数列の積を表示する
print product(5, 3, 10);
# 階乗を求める関数
sub Fact
{
my ($n) = @_;
if ($n <= 1)
{
1;
}
else
{
$n * Fact($n - 1);
}
}
# 10の階乗
print Fact(10), "\n";
print 10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1, "\n";
# 下降階乗冪
sub FallingFact
{
my ($x, $n) = @_;
if ($n <= 1)
{
$x;
}
else
{
$x * FallingFact($x - 1, $n - 1);
}
}
# 10 から 6 までの 総乗
print FallingFact(10, 5), "\n";
print 10 * 9 * 8 * 7 * 6, "\n";
# 上昇階乗冪
sub RisingFact
{
my ($x, $n) = @_;
if ($n <= 1)
{
$x;
}
else
{
$x * RisingFact($x + 1, $n - 1);
}
}
# 10 から 14 までの 総乗
print RisingFact(10, 5), "\n";
print 10 * 11 * 12 * 13 * 14, "\n";
# 階乗
sub Fact
{
my ($n) = @_;
if ($n <= 1)
{
1;
}
else
{
$n * Fact($n - 1);
}
}
# 下降階乗冪
sub FallingFact
{
my ($x, $n) = @_;
if ($n <= 1)
{
$x;
}
else
{
$x * FallingFact($x - 1, $n - 1);
}
}
# 順列 (異なる 10 個のものから 5 個取ってできる順列の総数)
$n = 10;
$r = 5;
print Fact($n) / Fact($n - $r), "\n";
print FallingFact($n, $r), "\n";
# 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数)
$n = 10;
$r = 5;
print $n ** $r, "\n";
# 組合せ
sub Comb
{
my ($n, $r) = @_;
if ($r == 0 || $r == $n)
{
1;
}
elsif ($r == 1)
{
$n;
}
else
{
Comb($n - 1, $r - 1) + Comb($n - 1, $r);
}
}
# 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数)
$n = 10;
$r = 5;
print Comb($n, $r), "\n";
# 組合せ
sub Comb
{
my ($n, $r) = @_;
if ($r == 0 || $r == $n)
{
1;
}
elsif ($r == 1)
{
$n;
}
else
{
Comb($n - 1, $r - 1) + Comb($n - 1, $r);
}
}
# 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数)
$n = 10;
$r = 5;
print Comb($n + $r - 1, $r), "\n";
use Math::Trig qw'deg2rad tan';
# 自作の正弦関数
sub mySin
{
my ($x, $n, $nega, $numerator, $denominator, $y) = @_;
$m = 2 * $n;
$denominator = $denominator * ($m + 1) * $m;
$numerator = $numerator * $x * $x;
$a = $numerator / $denominator;
# 十分な精度になったら処理を抜ける
if ($a <= 0.00000000001)
{
$y;
}
else
{
$y + mySin($x, ++$n, !$nega, $numerator, $denominator, $nega ? $a : -$a);
}
}
for $i (0..24)
{
$degree = $i * 15;
if ($degree % 30 == 0 || $degree % 45 == 0)
{
$radian = deg2rad($degree);
# 自作の正弦関数
$d1 = mySin($radian, 1, 0, $radian, 1.0, $radian);
# 標準の正弦関数
$d2 = sin($radian);
# 標準関数との差異
printf("%3d : %13.10f - %13.10f = %13.10f\n", $degree, $d1, $d2, $d1 - $d2);
}
}
use Math::Trig qw'deg2rad tan';
# 自作の余弦関数
sub myCos
{
my ($x, $n, $nega, $numerator, $denominator, $y) = @_;
$m = 2 * $n;
$denominator = $denominator * $m * ($m - 1);
$numerator = $numerator * $x * $x;
$a = $numerator / $denominator;
# 十分な精度になったら処理を抜ける
if ($a <= 0.00000000001)
{
$y;
}
else
{
$y + myCos($x, ++$n, !$nega, $numerator, $denominator, $nega ? $a : -$a);
}
}
for $i (0..24)
{
$degree = $i * 15;
if ($degree % 30 == 0 || $degree % 45 == 0)
{
$radian = deg2rad($degree);
# 自作の余弦関数
$d1 = myCos($radian, 1, 0, 1.0, 1.0, 1.0);
# 標準の余弦関数
$d2 = cos($radian);
# 標準関数との差異
printf("%3d : %13.10f - %13.10f = %13.10f\n", $degree, $d1, $d2, $d1 - $d2);
}
}
use Math::Trig qw'deg2rad tan';
# 自作の正接関数
sub myTan
{
my ($x, $x2, $n, $t) = @_;
$t = $x2 / ($n - $t);
$n -= 2;
if ($n <= 1)
{
$x / (1 - $t);
}
else
{
myTan($x, $x2, $n, $t);
}
}
for $i (0..12)
{
if (($i * 15) % 180 != 0)
{
$degree = $i * 15 - 90;
$radian = deg2rad($degree);
$x2 = $radian * $radian;
# 自作の正接関数
$d1 = myTan($radian, $x2, 15, 0.0); # 15:必要な精度が得られる十分大きな奇数
# 標準の正接関数
$d2 = tan($radian);
# 標準関数との差異
printf("%3d : %13.10f - %13.10f = %13.10f\n", $degree, $d1, $d2, $d1 - $d2);
}
}
# 自作の指数関数
sub myExp
{
my ($x, $n, $numerator, $denominator, $y) = @_;
$denominator = $denominator * $n;
$numerator = $numerator * $x;
$a = $numerator / $denominator;
# 十分な精度になったら処理を抜ける
if (abs($a) <= 0.00000000001)
{
$y;
}
else
{
$y + myExp($x, ++$n, $numerator, $denominator, $a);
}
}
for $i (0..20)
{
$x = ($i - 10) / 4.0;
# 標準の指数関数
$d1 = exp($x);
# 自作の指数関数
$d2 = myExp($x, 1, 1.0, 1.0, 1.0);
# 標準関数との差異
printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
# 自作の指数関数
sub myExp
{
my ($x, $x2, $n, $t) = @_;
$t = $x2 / ($n + $t);
$n -= 4;
if ($n < 6)
{
1 + ((2 * $x) / (2 - $x + $t));
}
else
{
myExp($x, $x2, $n, $t);
}
}
for $i (0..20)
{
$x = ($i - 10) / 4.0;
# 標準の指数関数
$d1 = exp($x);
# 自作の指数関数
$x2 = $x * $x;
$d2 = myExp($x, $x2, 30, 0.0); # 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる
# 標準関数との差異
printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
# 自作の対数関数
sub myLog
{
my ($x2, $numerator, $denominator, $y) = @_;
$denominator = $denominator + 2;
$numerator = $numerator * $x2 * $x2;
$a = $numerator / $denominator;
# 十分な精度になったら処理を抜ける
if (abs($a) <= 0.00000000001)
{
$y;
}
else
{
$y + myLog($x2, $numerator, $denominator, $a);
}
}
for $i (1..20)
{
$x = $i / 5.0;
# 標準の対数関数
$d1 = log($x);
# 自作の対数関数
$x2 = ($x - 1) / ($x + 1);
$d2 = 2 * myLog($x2, $x2, 1.0, $x2);
# 標準関数との差異
printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
# 自作の対数関数
sub myLog
{
my ($x, $n, $t) = @_;
$n2 = $n;
$x2 = $x;
if ($n > 3)
{
if ($n % 2 == 0)
{
$n2 = 2;
}
$x2 = $x * int($n / 2);
}
$t = $x2 / ($n2 + $t);
if ($n <= 2)
{
$x / (1 + $t);
}
else
{
myLog($x, --$n, $t);
}
}
for $i (1..20)
{
$x = $i / 5.0;
# 標準の対数関数
$d1 = log($x);
# 自作の対数関数
$d2 = myLog($x - 1, 27, 0.0); # 27:必要な精度が得られる十分大きな奇数
# 標準関数との差異
printf("%5.2f : %13.10f - %13.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
use Math::Trig;
# 自作の双曲線正弦関数
sub mySinh
{
my ($x, $n, $numerator, $denominator, $y) = @_;
$m = 2 * $n;
$denominator = $denominator * ($m + 1) * $m;
$numerator = $numerator * $x * $x;
$a = $numerator / $denominator;
# 十分な精度になったら処理を抜ける
if (abs($a) <= 0.00000000001)
{
$y;
}
else
{
$y + mySinh($x, ++$n, $numerator, $denominator, $a);
}
}
for $i (0..20)
{
$x = $i - 10;
# 自作の双曲線正弦関数
$d1 = mySinh($x, 1, $x, 1.0, $x);
# 標準の双曲線正弦関数
$d2 = sinh($x);
# 標準関数との差異
printf("%3d : %17.10f - %17.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
use Math::Trig;
# 自作の双曲線余弦関数
sub myCosh
{
my ($x, $n, $numerator, $denominator, $y) = @_;
$m = 2 * $n;
$denominator = $denominator * $m * ($m - 1);
$numerator = $numerator * $x * $x;
$a = $numerator / $denominator;
# 十分な精度になったら処理を抜ける
if (abs($a) <= 0.00000000001)
{
$y;
}
else
{
$y + myCosh($x, ++$n, $numerator, $denominator, $a);
}
}
for $i (0..20)
{
$x = $i - 10;
# 自作の双曲線余弦関数
$d1 = myCosh($x, 1, 1.0, 1.0, 1.0);
# 標準の双曲線余弦関数
$d2 = cosh($x);
# 標準関数との差異
printf("%3d : %17.10f - %17.10f = %13.10f\n", $x, $d1, $d2, $d1 - $d2);
}
use Math::Trig 'pi';
my $a = 0;
my $b = 1;
# 台形則で積分
my $n = 2;
for $j (1..10)
{
my $h = ($b - $a) / $n;
my $s = 0;
my $x = $a;
for $i (1..($n - 1))
{
$x += $h;
$s += f($x);
}
$s = $h * ((f($a) + f($b)) / 2 + $s);
$n *= 2;
# 結果を π と比較
printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - pi);
}
sub f
{
my ($x) = @_;
4 / (1 + $x * $x);
}
use Math::Trig 'pi';
my $a = 0;
my $b = 1;
# 中点則で積分
my $n = 2;
for $j (1..10)
{
my $h = ($b - $a) / $n;
my $s = 0;
my $x = $a + ($h / 2);
for $i (1..$n)
{
$s += f($x);
$x += $h;
}
$s *= $h;
$n *= 2;
# 結果を π と比較
printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - pi);
}
sub f
{
my ($x) = @_;
4 / (1 + $x * $x);
}
use Math::Trig 'pi';
my $a = 0;
my $b = 1;
# Simpson則で積分
my $n = 2;
for $j (1..5)
{
my $h = ($b - $a) / $n;
my $s2 = 0;
my $s4 = 0;
my $x = $a + $h;
for $i (1..($n / 2))
{
$s4 += f($x);
$x += $h;
$s2 += f($x);
$x += $h;
}
$s2 = ($s2 - f($b)) * 2 + f($a) + f($b);
$s4 *= 4;
my $s = ($s2 + $s4) * $h / 3;
$n *= 2;
# 結果を π と比較
printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - pi);
}
sub f
{
my ($x) = @_;
4 / (1 + $x * $x);
}
use Math::Trig 'pi';
my $a = 0;
my $b = 1;
my @t = ();
# 台形則で積分
my $n = 2;
for $i (1..6)
{
my $h = ($b - $a) / $n;
my $s = 0;
my $x = $a;
for $j (1..($n - 1))
{
$x += $h;
$s += f($x);
}
# 結果を保存
$t[$i][1] = $h * ((f($a) + f($b)) / 2 + $s);
$n *= 2;
}
# Richardsonの補外法
$n = 4;
for $j (2..6)
{
for $i ($j..6)
{
$t[$i][$j] = $t[$i][$j - 1] + ($t[$i][$j - 1] - $t[$i - 1][$j - 1]) / ($n - 1);
if ($i == $j)
{
# 結果を π と比較
printf("%2d : %13.10f, %13.10f\n", $j, $t[$i][$j], $t[$i][$j] - pi);
}
}
$n *= 4;
}
sub f
{
my ($x) = @_;
4 / (1 + $x * $x);
}
# データ点の数 - 1
use constant N => 6;
my @x = ();
my @y = ();
# 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for $i (0..N)
{
my $d = $i * 1.5 - 4.5;
$x[$i] = $d;
$y[$i] = f($d);
}
# 0.5刻みで 与えられていない値を補間
for $i (0..18)
{
my $d = $i * 0.5 - 4.5;
my $d1 = f($d);
my $d2 = lagrange($d, \@x, \@y);
# 元の関数と比較
printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d, $d1, $d2, $d1 - $d2);
}
# 元の関数
sub f
{
my ($x) = @_;
$x - ($x ** 3) / (3 * 2) + ($x ** 5) / (5 * 4 * 3 * 2);
}
# Lagrange (ラグランジュ) 補間
sub lagrange
{
my ($d, $x, $y) = @_;
my $sum = 0;
for $i (0..N)
{
my $prod = $$y[$i];
for $j (0..N)
{
if ($j != $i)
{
$prod *= ($d - $$x[$j]) / ($$x[$i] - $$x[$j]);
}
}
$sum += $prod;
}
$sum;
}
# データ点の数 - 1
use constant N => 6;
my @x = ();
my @y = ();
# 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for $i (0..N)
{
my $d = $i * 1.5 - 4.5;
$x[$i] = $d;
$y[$i] = f($d);
}
# 0.5刻みで 与えられていない値を補間
for $i (0..18)
{
my $d = $i * 0.5 - 4.5;
my $d1 = f($d);
my $d2 = neville($d, \@x, \@y);
# 元の関数と比較
printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d, $d1, $d2, $d1 - $d2);
}
# 元の関数
sub f
{
my ($x) = @_;
$x - ($x ** 3) / (3 * 2) + ($x ** 5) / (5 * 4 * 3 * 2);
}
# Neville (ネヴィル) 補間
sub neville
{
my ($d, $x, $y) = @_;
my @w = ();
for $i (0..N)
{
$w[0][$i] = $$y[$i];
}
for $j (1..N)
{
for $i (0..(N - $j))
{
$w[$j][$i] = $w[$j-1][$i+1] + ($w[$j-1][$i+1] - $w[$j-1][$i]) * ($d - $$x[$i+$j]) / ($$x[$i+$j] - $$x[$i]);
}
}
$w[N][0];
}
# データ点の数 - 1
use constant N => 6;
my @x = ();
my @y = ();
# 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for $i (0..N)
{
my $d1 = $i * 1.5 - 4.5;
$x[$i] = $d1;
$y[$i] = f($d1);
}
# 差分商の表を作る
my @d = ();
for $j (0..N)
{
$d[0][$j] = $y[$j];
}
for $i (1..N)
{
for $j (0..(N - $i))
{
$d[$i][$j] = ($d[$i-1][$j+1] - $d[$i-1][$j]) / ($x[$j+$i] - $x[$j]);
}
}
# ｎ階差分商
my @a = ();
for $j (0..N)
{
$a[$j] = $d[$j][0];
}
# 0.5刻みで 与えられていない値を補間
for $i (0..18)
{
my $d1 = $i * 0.5 - 4.5;
my $d2 = f($d1);
my $d3 = newton($d1, \@x, \@a);
# 元の関数と比較
printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d1, $d2, $d3, $d2 - $d3);
}
# 元の関数
sub f
{
my ($x) = @_;
$x - ($x ** 3) / (3 * 2) + ($x ** 5) / (5 * 4 * 3 * 2);
}
# Newton (ニュートン) 補間
sub newton
{
my ($d, $x, $a) = @_;
my $sum = $$a[0];
for $i (1..N)
{
my $prod = $$a[$i];
for $j (0..($i - 1))
{
$prod *= ($d - $$x[$j]);
}
$sum += $prod;
}
$sum;
}
# データ点の数 - 1
use constant N => 6;
use constant Nx2 => 13;
my @x = ();
my @y = ();
my @yd = ();
# 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for $i (0..N)
{
my $d1 = $i * 1.5 - 4.5;
$x[$i] = $d1;
$y[$i] = f($d1);
$yd[$i] = fd($d1);
}
# 差分商の表を作る
my @z = ();
my @d = ();
for $i (0..Nx2)
{
$j = int($i / 2);
$z[$i] = $x[$j];
$d[0][$i] = $y[$j];
}
for $i (1..Nx2)
{
for $j (0..(Nx2 - $i))
{
if ($i == 1 && $j % 2 == 0)
{
$d[$i][$j] = $yd[int($j / 2)];
}
else
{
$d[$i][$j] = ($d[$i-1][$j+1] - $d[$i-1][$j]) / ($z[$j+$i] - $z[$j]);
}
}
}
# ｎ階差分商
my @a = ();
for $j (0..Nx2)
{
$a[$j] = $d[$j][0];
}
# 0.5刻みで 与えられていない値を補間
for $i (0..18)
{
my $d1 = $i * 0.5 - 4.5;
my $d2 = f($d1);
my $d3 = hermite($d1, \@z, \@a);
# 元の関数と比較
printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d1, $d2, $d3, $d2 - $d3);
}
# 元の関数
sub f
{
my ($x) = @_;
$x - ($x ** 3) / (3 * 2) + ($x ** 5) / (5 * 4 * 3 * 2);
}
# 導関数
sub fd
{
my ($x) = @_;
1 - ($x ** 2) / 2 + ($x ** 4) / (4 * 3 * 2);
}
# Hermite (エルミート) 補間
sub hermite
{
my ($d, $z, $a) = @_;
my $sum = $$a[0];
for $i (1..Nx2)
{
my $prod = $$a[$i];
for $j (0..($i - 1))
{
$prod *= ($d - $$z[$j]);
}
$sum += $prod;
}
$sum;
}
# データ点の数 - 1
use constant N => 6;
my @x = ();
my @y = ();
# 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
for $i (0..N)
{
my $d1 = $i * 1.5 - 4.5;
$x[$i] = $d1;
$y[$i] = f($d1);
}
# ３項方程式の係数の表を作る
my @a = ();
my @b = ();
my @c = ();
my @d = ();
for $i (1..(N - 1))
{
$a[$i] = $x[$i] - $x[$i-1];
$b[$i] = 2.0 * ($x[$i+1] - $x[$i-1]);
$c[$i] = $x[$i+1] - $x[$i];
$d[$i] = 6.0 * (($y[$i+1] - $y[$i]) / ($x[$i+1] - $x[$i]) - ($y[$i] - $y[$i-1]) / ($x[$i] - $x[$i-1]));
}
# ３項方程式を解く (ト−マス法)
my @g = ();
my @s = ();
$g[1] = $b[1];
$s[1] = $d[1];
for $i (2..(N - 1))
{
$g[$i] = $b[$i] - $a[$i] * $c[$i-1] / $g[$i-1];
$s[$i] = $d[$i] - $a[$i] * $s[$i-1] / $g[$i-1];
}
my @z = ();
$z[0] = 0;
$z[N] = 0;
$z[N-1] = $s[N-1] / $g[N-1];
for ($i = N - 2; $i >= 1; $i--)
{
$z[$i] = ($s[$i] - $c[$i] * $z[$i+1]) / $g[$i];
}
# 0.5刻みで 与えられていない値を補間
for $i (0..18)
{
my $d1 = $i * 0.5 - 4.5;
my $d2 = f($d1);
my $d3 = spline($d1, \@x, \@y, \@z);
# 元の関数と比較
printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d1, $d2, $d3, $d2 - $d3);
}
# 元の関数
sub f
{
my ($x) = @_;
$x - ($x ** 3) / (3 * 2) + ($x ** 5) / (5 * 4 * 3 * 2);
}
# Spline (スプライン) 補間
sub spline
{
my ($d, $x, $y, $z) = @_;
# 補間関数値がどの区間にあるか
my $k = -1;
for $i (1..N)
{
if ($d <= $$x[$i])
{
$k = $i - 1;
last;
}
}
if ($k < 0)
{
$k = N;
}
$d1 = $$x[$k+1] - $d;
$d2 = $d - $$x[$k];
$d3 = $$x[$k+1] - $$x[$k];
($$z[$k] * ($d1 ** 3) + $$z[$k+1] * ($d2 ** 3)) / (6.0 * $d3) +
($$y[$k] / $d3 - $$z[$k] * $d3 / 6.0) * $d1 +
($$y[$k+1] / $d3 - $$z[$k+1] * $d3 / 6.0) * $d2;
}
use Math::Trig 'pi';
# 重力加速度
my $g = -9.8;
# 空気抵抗係数
my $k = -0.01;
# 時間間隔(秒)
my $h = 0.01;
# 角度
my $degree = 45;
my $radian = $degree * pi / 180.0;
# 初速 250 km/h -> 秒速に変換
my $v = int(250 * 1000 / 3600);
# 水平方向の速度
my @vx = ();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
my @vy = ();
$vy[0] = $v * sin($radian);
# 経過秒数
my $t = 0.0;
# 位置
my $x = 0.0;
my $y = 0.0;
# Euler法
for (my $i = 1; $y >= 0.0; $i++)
{
# 経過秒数
$t = $i * $h;
# 位置
$x += $h * $vx[0];
$y += $h * $vy[0];
printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", $t, $vx[0], $vy[0], $x, $y);
# 速度
$vx[1] = $vx[0] + $h * fx($vx[0], $vy[0]);
$vy[1] = $vy[0] + $h * fy($vx[0], $vy[0]);
$vx[0] = $vx[1];
$vy[0] = $vy[1];
}
# 空気抵抗による水平方向の減速分
sub fx
{
my ($vx, $vy) = @_;
$k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
sub fy
{
my ($vx, $vy) = @_;
$g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
use Math::Trig 'pi';
# 重力加速度
my $g = -9.8;
# 空気抵抗係数
my $k = -0.01;
# 時間間隔(秒)
my $h = 0.01;
# 角度
my $degree = 45;
my $radian = $degree * pi / 180.0;
# 初速 250 km/h -> 秒速に変換
my $v = int(250 * 1000 / 3600);
# 水平方向の速度
my @vx = ();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
my @vy = ();
$vy[0] = $v * sin($radian);
# 経過秒数
my $t = 0.0;
# 位置
my @x = ();
$x[0] = 0.0;
my @y = ();
$y[0] = 0.0;
# Heun法
for (my $i = 1; $y[0] >= 0.0; $i++)
{
# 経過秒数
$t = $i * $h;
# 位置・速度
$x[1] = $x[0] + $h * $vx[0];
$y[1] = $y[0] + $h * $vy[0];
$vx[1] = $vx[0] + $h * fx($vx[0], $vy[0]);
$vy[1] = $vy[0] + $h * fy($vx[0], $vy[0]);
$x[2] = $x[0] + $h * ( $vx[0] + $vx[1] ) / 2;
$y[2] = $y[0] + $h * ( $vy[0] + $vy[1] ) / 2;
$vx[2] = $vx[0] + $h * (fx($vx[0], $vy[0]) + fx($vx[1], $vy[1])) / 2;
$vy[2] = $vy[0] + $h * (fy($vx[0], $vy[0]) + fy($vx[1], $vy[1])) / 2;
$x[0] = $x[2];
$y[0] = $y[2];
$vx[0] = $vx[2];
$vy[0] = $vy[2];
printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}
# 空気抵抗による水平方向の減速分
sub fx
{
my ($vx, $vy) = @_;
$k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
sub fy
{
my ($vx, $vy) = @_;
$g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
use Math::Trig 'pi';
# 重力加速度
my $g = -9.8;
# 空気抵抗係数
my $k = -0.01;
# 時間間隔(秒)
my $h = 0.01;
# 角度
my $degree = 45;
my $radian = $degree * pi / 180.0;
# 初速 250 km/h -> 秒速に変換
my $v = int(250 * 1000 / 3600);
# 水平方向の速度
my @vx = ();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
my @vy = ();
$vy[0] = $v * sin($radian);
# 経過秒数
my $t = 0.0;
# 位置
my @x = ();
$x[0] = 0.0;
my @y = ();
$y[0] = 0.0;
# 中点法
for (my $i = 1; $y[0] >= 0.0; $i++)
{
# 経過秒数
$t = $i * $h;
# 位置・速度
$vx[1] = $h * fx($vx[0], $vy[0]);
$vy[1] = $h * fy($vx[0], $vy[0]);
my $wx = $vx[0] + $vx[1] / 2;
my $wy = $vy[0] + $vy[1] / 2;
$vx[0] = $vx[0] + $h * fx($wx, $wy);
$vy[0] = $vy[0] + $h * fy($wx, $wy);
$x[0] = $x[0] + $h * $wx;
$y[0] = $y[0] + $h * $wy;
printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}
# 空気抵抗による水平方向の減速分
sub fx
{
my ($vx, $vy) = @_;
$k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
sub fy
{
my ($vx, $vy) = @_;
$g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
use Math::Trig 'pi';
# 重力加速度
my $g = -9.8;
# 空気抵抗係数
my $k = -0.01;
# 時間間隔(秒)
my $h = 0.01;
# 角度
my $degree = 45;
my $radian = $degree * pi / 180.0;
# 初速 250 km/h -> 秒速に変換
my $v = int(250 * 1000 / 3600);
# 水平方向の速度
my @vx = ();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
my @vy = ();
$vy[0] = $v * sin($radian);
# 経過秒数
my $t = 0.0;
# 位置
my @x = ();
$x[0] = 0.0;
my @y = ();
$y[0] = 0.0;
# Runge-Kutta法
for (my $i = 1; $y[0] >= 0.0; $i++)
{
# 経過秒数
$t = $i * $h;
# 位置・速度
$x[1] = $h * $vx[0];
$y[1] = $h * $vy[0];
$vx[1] = $h * fx($vx[0], $vy[0]);
$vy[1] = $h * fy($vx[0], $vy[0]);
my $wx = $vx[0] + $vx[1] / 2;
my $wy = $vy[0] + $vy[1] / 2;
$x[2] = $h * $wx;
$y[2] = $h * $wy;
$vx[2] = $h * fx($wx, $wy);
$vy[2] = $h * fy($wx, $wy);
$wx = $vx[0] + $vx[2] / 2;
$wy = $vy[0] + $vy[2] / 2;
$x[3] = $h * $wx;
$y[3] = $h * $wy;
$vx[3] = $h * fx($wx, $wy);
$vy[3] = $h * fy($wx, $wy);
$wx = $vx[0] + $vx[3];
$wy = $vy[0] + $vy[3];
$x[4] = $h * $wx;
$y[4] = $h * $wy;
$vx[4] = $h * fx($wx, $wy);
$vy[4] = $h * fy($wx, $wy);
$x[0] += ( $x[1] + $x[2] * 2 + $x[3] * 2 + $x[4]) / 6;
$y[0] += ( $y[1] + $y[2] * 2 + $y[3] * 2 + $y[4]) / 6;
$vx[0] += ($vx[1] + $vx[2] * 2 + $vx[3] * 2 + $vx[4]) / 6;
$vy[0] += ($vy[1] + $vy[2] * 2 + $vy[3] * 2 + $vy[4]) / 6;
printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}
# 空気抵抗による水平方向の減速分
sub fx
{
my ($vx, $vy) = @_;
$k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
sub fy
{
my ($vx, $vy) = @_;
$g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
use Math::Trig 'pi';
# 重力加速度
my $g = -9.8;
# 空気抵抗係数
my $k = -0.01;
# 時間間隔(秒)
my $h = 0.01;
# 角度
my $degree = 45;
my $radian = $degree * pi / 180.0;
# 初速 250 km/h -> 秒速に変換
my $v = int(250 * 1000 / 3600);
# 水平方向の速度
my @vx = ();
$vx[0] = $v * cos($radian);
# 鉛直方向の速度
my @vy = ();
$vy[0] = $v * sin($radian);
# 経過秒数
my $t = 0.0;
# 位置
my @x = ();
$x[0] = 0.0;
my @y = ();
$y[0] = 0.0;
# Runge-Kutta-Gill法
for (my $i = 1; $y[0] >= 0.0; $i++)
{
# 経過秒数
$t = $i * $h;
# 位置・速度
$x[1] = $h * $vx[0];
$y[1] = $h * $vy[0];
$vx[1] = $h * fx($vx[0], $vy[0]);
$vy[1] = $h * fy($vx[0], $vy[0]);
my $wx = $vx[0] + $vx[1] / 2;
my $wy = $vy[0] + $vy[1] / 2;
$x[2] = $h * $wx;
$y[2] = $h * $wy;
$vx[2] = $h * fx($wx, $wy);
$vy[2] = $h * fy($wx, $wy);
$wx = $vx[0] + $vx[1] * ((sqrt(2.0) - 1) / 2) + $vx[2] * (1 - 1 / sqrt(2.0));
$wy = $vy[0] + $vy[1] * ((sqrt(2.0) - 1) / 2) + $vy[2] * (1 - 1 / sqrt(2.0));
$x[3] = $h * $wx;
$y[3] = $h * $wy;
$vx[3] = $h * fx($wx, $wy);
$vy[3] = $h * fy($wx, $wy);
$wx = $vx[0] - $vx[2] / sqrt(2.0) + $vx[3] * (1 + 1 / sqrt(2.0));
$wy = $vy[0] - $vy[2] / sqrt(2.0) + $vy[3] * (1 + 1 / sqrt(2.0));
$x[4] = $h * $wx;
$y[4] = $h * $wy;
$vx[4] = $h * fx($wx, $wy);
$vy[4] = $h * fy($wx, $wy);
$x[0] += ( $x[1] + $x[2] * (2 - sqrt(2.0)) + $x[3] * (2 + sqrt(2.0)) + $x[4]) / 6;
$y[0] += ( $y[1] + $y[2] * (2 - sqrt(2.0)) + $y[3] * (2 + sqrt(2.0)) + $y[4]) / 6;
$vx[0] += ($vx[1] + $vx[2] * (2 - sqrt(2.0)) + $vx[3] * (2 + sqrt(2.0)) + $vx[4]) / 6;
$vy[0] += ($vy[1] + $vy[2] * (2 - sqrt(2.0)) + $vy[3] * (2 + sqrt(2.0)) + $vy[4]) / 6;
printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f\n", $t, $vx[0], $vy[0], $x[0], $y[0]);
}
# 空気抵抗による水平方向の減速分
sub fx
{
my ($vx, $vy) = @_;
$k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
sub fy
{
my ($vx, $vy) = @_;
$g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}
my $a = 1;
my $b = 2;
printf("%12.10f\n", bisection($a, $b));
sub bisection
{
my ($a, $b) = @_;
my $c;
while (1)
{
# 区間 (a, b) の中点 c = (a + b) / 2
$c = ($a + $b) / 2;
printf("%12.10f\t%13.10f\n", $c, $c - sqrt(2));
my $fc = f($c);
if (abs($fc) < 0.0000000001)
{
last;
}
if ($fc < 0)
{
# f(c) < 0 であれば, 解は区間 (c, b) の中に存在
$a = $c;
}
else
{
# f(c) > 0 であれば, 解は区間 (a, c) の中に存在
$b = $c;
}
}
$c;
}
sub f
{
my ($x) = @_;
$x * $x - 2;
}
my $a = 1;
my $b = 2;
printf("%12.10f\n", falseposition($a, $b));
sub falseposition
{
my ($a, $b) = @_;
my $c;
while (1)
{
# 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点
$c = ($a * f($b) - $b * f($a)) / (f($b) - f($a));
printf("%12.10f\t%13.10f\n", $c, $c - sqrt(2));
my $fc = f($c);
if (abs($fc) < 0.0000000001)
{
last;
}
if ($fc < 0)
{
# f(c) < 0 であれば, 解は区間 (c, b) の中に存在
$a = $c;
}
else
{
# f(c) > 0 であれば, 解は区間 (a, c) の中に存在
$b = $c;
}
}
$c;
}
sub f
{
my ($x) = @_;
$x * $x - 2;
}
my $x = 1;
printf("%12.10f\n", iterative($x));
sub iterative
{
my ($x0) = @_;
my $x1;
while (1)
{
$x1 = g($x0);
printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));
if (abs($x1 - $x0) < 0.0000000001)
{
last;
}
$x0 = $x1;
}
$x1;
}
sub g
{
my ($x) = @_;
($x / 2) + (1 / $x);
}
my $x = 2;
printf("%12.10f\n", newton($x));
sub newton
{
my ($x0) = @_;
my $x1;
while (1)
{
$x1 = $x0 - (f0($x0) / f1($x0));
printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));
if (abs($x1 - $x0) < 0.0000000001)
{
last;
}
$x0 = $x1;
}
$x1;
}
sub f0
{
my ($x) = @_;
$x * $x - 2;
}
sub f1
{
my ($x) = @_;
2 * $x;
}
my $x = 2;
printf("%12.10f\n", bailey($x));
sub bailey
{
my ($x0) = @_;
my $x1;
while (1)
{
$x1 = $x0 - (f0($x0) / (f1($x0) - (f0($x0) * f2($x0) / (2 * f1($x0)))));
printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));
if (abs($x1 - $x0) < 0.0000000001)
{
last;
}
$x0 = $x1;
}
$x1;
}
sub f0
{
my ($x) = @_;
$x * $x - 2;
}
sub f1
{
my ($x) = @_;
2 * $x;
}
sub f2
{
my ($x) = @_;
2;
}
my $x0 = 1;
my $x1 = 2;
printf("%12.10f\n", secant($x0, $x1));
sub secant
{
my ($x0, $x1) = @_;
my $x2;
while (1)
{
$x2 = $x1 - f($x1) * ($x1 - $x0) / (f($x1) - f($x0));
printf("%12.10f\t%13.10f\n", $x2, $x2 - sqrt(2));
if (abs($x2 - $x1) < 0.0000000001)
{
last;
}
$x0 = $x1;
$x1 = $x2;
}
$x2;
}
sub f
{
my ($x) = @_;
$x * $x - 2;
}
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];
if (abs($x1[$i] - $$x0[$i]) > 0.0000000001) {$finish = 0;}
}
for $i (0..N)
{
$$x0[$i] = $x1[$i];
printf("%14.10f\t", $$x0[$i]);
}
print "\n";
last 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];
if (abs($x1 - $$x0[$i]) > 0.0000000001) {$finish = 0;}
$$x0[$i] = $x1;
}
for $i (0..N)
{
printf("%14.10f\t", $$x0[$i]);
}
print "\n";
last 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";
disp_progress(\@a,\@b);
# 前進消去
forward_elimination(\@a,\@b);
print "前進消去後\n";
disp_progress(\@a,\@b);
# 後退代入
backward_substitution(\@a,\@b);
print "解\n";
for $i (0..N)
{
printf("%14.10f\t", $b[$i]);
}
print "\n";
# 前進消去
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_progress
{
my ($a, $b) = @_;
for $i (0..N)
{
for $j (0..N)
{
printf("%14.10f\t", $a[$i][$j]);
}
printf("%14.10f\n", $b[$i]);
}
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";
disp_progress(\@a,\@b);
# 前進消去
forward_elimination(\@a,\@b);
print "前進消去後\n";
disp_progress(\@a,\@b);
# 後退代入
backward_substitution(\@a,\@b);
print "解\n";
for $i (0..N)
{
printf("%14.10f\t", $b[$i]);
}
print "\n";
# 前進消去
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_progress
{
my ($a, $b) = @_;
for $i (0..N)
{
for $j (0..N)
{
printf("%14.10f\t", $a[$i][$j]);
}
printf("%14.10f\n", $b[$i]);
}
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";
disp_progress(\@a,\@b);
# 前進消去
forward_elimination(\@a,\@b);
print "前進消去後\n";
disp_progress(\@a,\@b);
# 後退代入
backward_substitution(\@a,\@b);
print "解\n";
for $i (0..N)
{
printf("%14.10f\t", $b[$i]);
}
print "\n";
# 前進消去
sub forward_elimination
{
my ($a, $b) = @_;
# 対角線上の係数を 1 にする
for $pivot (0..N)
{
my $s = $$a[$pivot][$pivot];
for $col (0..N)
{
$$a[$pivot][$col] /= $s;
}
$$b[$pivot] /= $s;
}
print "対角線上の係数を 1 にする\n";
disp_progress($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_progress
{
my ($a, $b) = @_;
for $i (0..N)
{
for $j (0..N)
{
printf("%14.10f\t", $a[$i][$j]);
}
printf("%14.10f\n", $b[$i]);
}
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";
disp_progress(\@a,\@b);
# LU分解
my @x = (0,0,0,0);
decomp(\@a,\@b,\@x);
print "X\n";
for $i (0..N)
{
printf("%14.10f\t", $x[$i]);
}
print "\n";
# LU分解
sub decomp
{
my ($a, $b, $x) = @_;
# 前進消去
forward_elimination(\@a,\@b);
# 下三角行列を確認
print "L\n";
for $i (0..N)
{
for $j (0..N)
{
if ($i > $j)
{
printf("%14.10f\t", $a[$i][$j]);
}
elsif ($i == $j)
{
printf("%14.10f\t", 1.0);
}
else
{
printf("%14.10f\t", 0.0);
}
}
print "\n";
}
print "\n";
# 上三角行列を確認
print "U\n";
for $i (0..N)
{
for $j (0..N)
{
if ($i <= $j)
{
printf("%14.10f\t", $a[$i][$j]);
}
else
{
printf("%14.10f\t", 0.0);
}
}
print "\n";
}
print "\n";
# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);
# y の 値 を確認
print "Y\n";
for $i (0..N)
{
printf("%14.10f\t", $y[$i]);
}
print "\n";
# Ux=y から x を求める (後退代入)
backward_substitution(\@a,\@y,\@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_progress
{
my ($a, $b) = @_;
for $i (0..N)
{
for $j (0..N)
{
printf("%14.10f\t", $a[$i][$j]);
}
printf("%14.10f\n", $b[$i]);
}
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_progress(\@a);
# LL^T分解
my @x = (0,0,0,0);
decomp(\@a,\@b,\@x);
print "X\n";
for $i (0..N)
{
printf("%14.10f\t", $x[$i]);
}
print "\n";
# LL^T分解
sub decomp
{
my ($a, $b, $x) = @_;
# 前進消去
forward_elimination(\@a,\@b);
print "L と L^T\n";
for $i (0..N)
{
for $j (0..N)
{
if ($j <= $i)
{
printf("%14.10f\t", $a[$i][$j]);
}
else
{
printf("%14.10f\t", $a[$j][$i]);
}
}
print "\n";
}
print "\n";
# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);
# y の 値 を確認
print "Y\n";
for $i (0..N)
{
printf("%14.10f\t", $y[$i]);
}
print "\n";
# Ux=y から x を求める (後退代入)
backward_substitution(\@a,\@y,\@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];
}
}
}
# 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[$col][$row] * $$x[$col];
}
$$x[$row] = $$y[$row] / $a[$row][$row];
}
}
# 状態を確認
sub disp_progress
{
my ($a) = @_;
for $i (0..N)
{
for $j (0..N)
{
printf("%14.10f\t", $a[$i][$j]);
}
print "\n";
}
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_progress(\@a);
# LDL^T分解
my @x = (0,0,0,0);
decomp(\@a,\@b,\@x);
print "X\n";
for $i (0..N)
{
printf("%14.10f\t", $x[$i]);
}
print "\n";
# LDL^T分解
sub decomp
{
my ($a, $b, $x) = @_;
# 前進消去
forward_elimination(\@a,\@b);
print "L と D\n";
for $i (0..N)
{
for $j (0..N)
{
if ($j <= $i)
{
printf("%14.10f\t", $a[$i][$j]);
}
else
{
printf("%14.10f\t", $a[$j][$i]);
}
}
print "\n";
}
print "\n";
# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);
# y の 値 を確認
print "Y\n";
for $i (0..N)
{
printf("%14.10f\t", $y[$i]);
}
print "\n";
# DL^Tx=y から x を求める (後退代入)
backward_substitution(\@a,\@y,\@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];
}
# 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[$col][$row] * $a[$row][$row] * $$x[$col];
}
$$x[$row] = $$y[$row] / $a[$row][$row];
}
}
# 状態を確認
sub disp_progress
{
my ($a) = @_;
for $i (0..N)
{
for $j (0..N)
{
printf("%14.10f\t", $a[$i][$j]);
}
print "\n";
}
print "\n";
}
|