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 までの 自然数の2乗の和 my $n = 1000; print($n * ($n + 1) * (2 * $n + 1) / 6, "\n");
# 100 までの 自然数の3乗の和 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 まで, 7点だけ値をセット 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 まで, 7点だけ値をセット 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 まで, 7点だけ値をセット 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]); } } # n階差分商 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 まで, 7点だけ値をセット 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]); } } } # n階差分商 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 まで, 7点だけ値をセット for $i (0..N) { my $d1 = $i * 1.5 - 4.5; $x[$i] = $d1; $y[$i] = f($d1); } # 3項方程式の係数の表を作る 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])); } # 3項方程式を解く (ト−マス法) 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 "X\n"; disp_vector(\@c); # ヤコビの反復法 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]; } return if ($finish); disp_vector($x0); } } # 1次元配列を表示 sub disp_vector { my ($row) = @_; foreach $col (@$row) { printf("%14.10f\t", $col); } print "\n"; }
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 "X\n"; disp_vector(\@c); # ガウス・ザイデル法 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; } return if ($finish); disp_vector($x0); } } # 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 "pivoting\n"; print "A\n"; disp_matrix(\@a); print "B\n"; disp_vector(\@b); print "\n"; # 前進消去 forward_elimination(\@a,\@b); print "forward elimination\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 "pivoting\n"; print "A\n"; disp_matrix(\@a); print "B\n"; disp_vector(\@b); print "\n"; # 前進消去 forward_elimination(\@a,\@b); print "forward elimination\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 "pivoting\n"; print "A\n"; disp_matrix(\@a); print "B\n"; disp_vector(\@b); print "\n"; # 前進消去 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"; # 前進消去 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); # L^Tx=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]; } } } # 前進代入 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]; } } # 後退代入 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"; # 前進消去 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); # DL^Tx=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; } } # 前進代入 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]; } } # 後退代入 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"; }
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 @x = (1.0, 0.0, 0.0, 0.0); # ベキ乗法 Power iteration $lambda = power(\@a, \@x); print "\neigenvalue\n"; printf("%14.10f", $lambda); print "\neigenvector\n"; disp_vector(\@x); } # ベキ乗法 Power iteration sub power { my ($a, $x0) = @_; my $lambda = 0.0; # 正規化 (ベクトル x0 の長さを1にする) normarize(\@$x0); my $e0 = 0.0; for $i (0..N) { $e0 += $$x0[$i]; } for $k (1..200) { # 1次元配列を表示 printf("%3d\t", $k); disp_vector(\@$x0); # 行列の積 x1 = A × x0 my @x1; for $i (0..N) { for $j (0..N) { $x1[$i] += $$a[$i][$j] * $$x0[$j]; } } # 内積 inner product my $p0 = 0.0; my $p1 = 0.0; for $i (0..N) { $p0 += $x1[$i] * $x1[$i]; $p1 += $x1[$i] * $$x0[$i]; } # 固有値 eigenvalue $lambda = $p0 / $p1; # 正規化 (ベクトル x1 の長さを1にする) normarize(\@x1); # 収束判定 my $e1 = 0.0; for $i (0..N) { $e1 += $x1[$i]; } last if (abs($e0 - $e1) < 0.00000000001); for $i (0..N) { $$x0[$i] = $x1[$i]; } $e0 = $e1; } $lambda; } # 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; } }
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 @x = (1.0, 0.0, 0.0, 0.0); # LU分解 forward_elimination(\@a); # 逆ベキ乗法 Inverse iteration $lambda = inverse(\@a, \@x); print "\neigenvalue\n"; printf("%14.10f", $lambda); print "\neigenvector\n"; disp_vector(\@x); } # 逆ベキ乗法 Inverse iteration sub inverse { my ($a, $x0) = @_; my $lambda = 0.0; # 正規化 (ベクトル x0 の長さを1にする) normarize(\@$x0); my $e0 = 0.0; for $i (0..N) { $e0 += $$x0[$i]; } for $k (1..200) { # 1次元配列を表示 printf("%3d\t", $k); disp_vector(\@$x0); # Ly = b から y を求める (前進代入) my @b = (0.0, 0.0, 0.0, 0.0); my @y; for $i (0..N) { $b[$i] = $$x0[$i]; } forward_substitution(\@$a, \@y, \@b); # Ux = y から x を求める (後退代入) my @x1; backward_substitution(\@$a, \@x1, \@y); # 内積 inner product my $p0 = 0.0; my $p1 = 0.0; for $i (0..N) { $p0 += $x1[$i] * $x1[$i]; $p1 += $x1[$i] * $$x0[$i]; } # 固有値 eigenvalue $lambda = $p1 / $p0; # 正規化 (ベクトル x1 の長さを1にする) normarize(\@x1); # 収束判定 my $e1 = 0.0; for $i (0..N) { $e1 += $x1[$i]; } last if (abs($e0 - $e1) < 0.00000000001); for $i (0..N) { $$x0[$i] = $x1[$i]; } $e0 = $e1; } $lambda; } # LU分解 sub forward_elimination { my ($a) = @_; 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; # これが 下三角行列 } } } # 前進代入 sub forward_substitution { my ($a, $y, $b) = @_; for $row (0..N) { for $col (0..$row) { $$b[$row] -= $$a[$row][$col] * $$y[$col]; } $$y[$row] = $$b[$row]; } } # 後退代入 sub backward_substitution { my ($a, $x, $y) = @_; for $row (reverse(0..N)) { for $col (reverse(($row + 1)..N)) { $$y[$row] -= $$a[$row][$col] * $$x[$col]; } $$x[$row] = $$y[$row] / $$a[$row][$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; } }
use constant N => 3; # LR分解で固有値を求める 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 @l, @u; for $k (1..200) { # LU分解 decomp(\@a, \@l, \@u); # 行列の積 multiply(\@u, \@l, \@a); # 対角要素を表示 printf("%3d\t", $k); disp_eigenvalue(\@a); # 収束判定 my $e = 0.0; for $i (1..N) { for $j (0..($i - 1)) { $e += abs($a[$i][$j]); } } last if ($e < 0.00000000001); } print "\neigenvalue\n"; disp_eigenvalue(\@a); } # LU分解 sub decomp { my ($a, $l, $u) = @_; for $i (0..N) { for $j (1..N) { $$l[$i][$j] = 0.0; $$u[$i][$j] = 0.0; } } $$l[0][0] = 1.0; for $j (0..N) { $$u[0][$j] = $$a[0][$j]; } for $i (1..N) { $$u[$i][0] = 0.0; $$l[0][$i] = 0.0; $$l[$i][0] = $$a[$i][0] / $$u[0][0]; } for $i (1..N) { $$l[$i][$i] = 1.0; $t = $$a[$i][$i]; for $k (0..$i) { $t -= $$l[$i][$k] * $$u[$k][$i]; } $$u[$i][$i] = $t; if ($i < N) { for $j (($i + 1)..N) { $$u[$j][$i] = 0; $$l[$i][$j] = 0; $t = $$a[$j][$i]; for $k (0..$i) { $t -= $$l[$j][$k] * $$u[$k][$i]; } $$l[$j][$i] = $t / $$u[$i][$i]; $t = $$a[$i][$j]; for $k (0..$i) { $t -= $$l[$i][$k] * $$u[$k][$j]; } $$u[$i][$j] = $t; } } } } # 行列の積 sub multiply { my ($a, $b, $c) = @_; for $i (0..N) { for $j (0..N) { $s = 0.0; for $k (0..N) { $s += $$a[$i][$k] * $$b[$k][$j]; } $$c[$i][$j] = $s; } } } # 対角要素を表示 sub disp_eigenvalue { my ($a) = @_; for $i (0..N) { printf("%14.10f\t", $$a[$i][$i]); } print("\n"); }
use constant N => 3; # QR分解で固有値を求める main(); sub main { my @a = ([5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]); my @q, @r; for $k (1..100) { # QR分解 decomp(\@a, \@q, \@r); # 行列の積 multiply(\@r, \@q, \@a); # 対角要素を表示 printf("%3d\t", $k); disp_eigenvalue(\@a); # 収束判定 my $e = 0.0; for $i (1..N) { for $j (0..($i - 1)) { $e += abs($a[$i][$j]); } } last if ($e < 0.00000000001); } print "\neigenvalue\n"; disp_eigenvalue(\@a); } # QR分解 sub decomp { my ($a, $q, $r) = @_; my @x; for $k (0..N) { for $i (0..N) { $x[$i] = $$a[$i][$k]; } for $j (0..($k - 1)) { my $t = 0.0; for $i (0..N) { $t += $$a[$i][$k] * $$q[$i][$j]; } $$r[$j][$k] = $t; $$r[$k][$j] = 0.0; for $i (0..N) { $x[$i] -= $t * $$q[$i][$j]; } } my $s = 0.0; for $i (0..N) { $s += $x[$i] * $x[$i]; } $$r[$k][$k] = sqrt($s); for $i (0..N) { $$q[$i][$k] = $x[$i] / $$r[$k][$k]; } } } # 行列の積 sub multiply { my ($a, $b, $c) = @_; for $i (0..N) { for $j (0..N) { $s = 0.0; for $k (0..N) { $s += $$a[$i][$k] * $$b[$k][$j]; } $$c[$i][$j] = $s; } } } # 対角要素を表示 sub disp_eigenvalue { my ($a) = @_; for $i (0..N) { printf("%14.10f\t", $$a[$i][$i]); } print("\n"); }
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; } }
use constant N => 3; # ハウスホルダー変換とQR分解で固有値を求める main(); sub main { my @a = ([5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]); my @d = (0.0, 0.0, 0.0, 0.0); my @e = (0.0, 0.0, 0.0, 0.0); # ハウスホルダー変換 tridiagonalize(\@a, \@d, \@e); # QR分解 decomp(\@a, \@d, \@e); print "\neigenvalue\n"; disp_vector(\@d); print "\neigenvector\n"; disp_matrix(\@a); } # ハウスホルダー変換 sub tridiagonalize { my ($a, $d, $e) = @_; my @v = (0.0, 0.0, 0.0, 0.0); for $k (0..(N - 2)) { my @w = (0.0, 0.0, 0.0, 0.0); $$d[$k] = $$a[$k][$k]; my $t = 0.0; for $i (($k + 1)..N) { $w[$i] = $$a[$i][$k]; $t += ($w[$i] * $w[$i]); } my $s = sqrt($t); $s = -$s if ($w[$k + 1] < 0); if (abs($s) < 0.00000000001) { $$e[$k + 1] = 0.0; } else { $$e[$k + 1] = -$s; $w[$k + 1] += $s; $s = 1 / sqrt($w[$k + 1] * $s); for $i (($k + 1)..N) { $w[$i] *= $s; } for $i (($k + 1)..N) { $s = 0.0; for $j (($k + 1)..N) { if ($j <= $i) { $s += $$a[$i][$j] * $w[$j]; } else { $s += $$a[$j][$i] * $w[$j]; } } $v[$i] = $s; } $s = 0; for $i (($k + 1)..N) { $s += $w[$i] * $v[$i]; } $s /= 2; for $i (($k + 1)..N) { $v[$i] -= $s * $w[$i]; } for $i (($k + 1)..N) { for $j (($k + 1)..$i) { $$a[$i][$j] -= $w[$i] * $v[$j] + $w[$j] * $v[$i]; } } # 次の行は固有ベクトルを求めないなら不要 for $i (($k + 1)..N) { $$a[$i][$k] = $w[$i]; } } } $$d[N - 1] = $$a[N - 1][N - 1]; $$d[N] = $$a[N][N]; $$e[0] = 0.0; $$e[N] = $$a[N][N - 1]; # 次の行は固有ベクトルを求めないなら不要 for $k (reverse(0..N)) { $w = (0.0, 0.0, 0.0, 0.0); if ($k < (N - 1)) { for $i (($k + 1)..N) { $w[$i] = $$a[$i][$k]; } for $i (($k + 1)..N) { $s = 0.0; for $j (($k + 1)..N) { $s += $$a[$i][$j] * $w[$j]; } $v[$i] = $s; } for $i (($k + 1)..N) { for $j (($k + 1)..N) { $$a[$i][$j] -= $v[$i] * $w[$j]; } } } for $i (0..N) { $$a[$i][$k] = 0.0; } $$a[$k][$k] = 1.0; } } # QR分解 sub decomp { my ($a, $d, $e) = @_; $$e[0] = 1.0; my $h = N; $h-- while (abs($$e[$h]) < 0.00000000001); while ($h > 0) { $$e[0] = 0.0; my $l = $h - 1; $l-- while (abs($$e[$l]) >= 0.00000000001); for $j (1..100) { my $w = ($$d[$h - 1] - $$d[$h]) / 2.0; my $s = sqrt($w * $w + $$e[$h] * $$e[$h]); $s = - $s if ($w < 0.0); my $x = $$d[$l] - $$d[$h] + $$e[$h] * $$e[$h] / ($w + $s); my $y = $$e[$l + 1]; my $z = 0.0; my $t = 0.0; my $u = 0.0; for $k ($l..($h - 1)) { if (abs($x) >= abs($y)) { $t = -$y / $x; $u = 1 / sqrt($t * $t + 1.0); $s = $t * $u; } else { $t = -$x / $y; $s = 1 / sqrt($t * $t + 1.0); $s = -$s if ($t < 0.0); $u = $t * $s; } $w = $$d[$k] - $$d[$k + 1]; $t = ($w * $s + 2 * $u * $$e[$k + 1]) * $s; $$d[$k ] = $$d[$k ] - $t; $$d[$k + 1] = $$d[$k + 1] + $t; $$e[$k ] = $u * $$e[$k] - $s * $z; $$e[$k + 1] = $$e[$k + 1] * ($u * $u - $s * $s) + $w * $s * $u; # 次の行は固有ベクトルを求めないなら不要 for $i (0..N) { $x = $$a[$k ][$i]; $y = $$a[$k + 1][$i]; $$a[$k ][$i] = $u * $x - $s * $y; $$a[$k + 1][$i] = $s * $x + $u * $y; } if ($k < N - 1) { $x = $$e[$k + 1]; $y = -$s * $$e[$k + 2]; $z = $y; $$e[$k + 2] = $u * $$e[$k + 2]; } } printf("%3d\t", $j); disp_vector(\@$d); # 収束判定 last if (abs($$e[$h]) < 0.00000000001); } $$e[0] = 1.0; $h-- while (abs($$e[$h]) < 0.00000000001); } # 次の行は固有ベクトルを求めないなら不要 for $k (0..(N - 1)) { $l = $k; for $i (($k + 1)..N) { $l = $i if ($$d[$i] > $$d[$l]); } $t = $$d[$k]; $$d[$k] = $$d[$l]; $$d[$l] = $t; for $i (0..N) { $t = $$a[$k][$i]; $$a[$k][$i] = $$a[$l][$i]; $$a[$l][$i] = $t; } } } # 1次元配列を表示 sub disp_vector { my ($x) = @_; for $i (0..N) { printf("%14.10f\t", $$x[$i]); } print "\n"; } # 2次元配列を表示 sub disp_matrix { my ($matrix) = @_; foreach $row (@$matrix) { foreach $col (@$row) { printf("%14.10f\t", $col); } print "\n"; } }