Echo (3 + 5) Echo (3 - 5) Echo (3 * 5) Echo ([Math]::Pow(3, 5)) Echo (5 / 3) Echo ([Math]::Floor(5 / 3)) Echo (5 % 3) Echo (5 % 3) Write-Output (3 * 5) Write-Host (3 * 5) Write-Host ([String]::Format("{0,3:D}`n", 3 * 5)) -nonewline Write-Host ([String]::Format("{0,23:F20}", 5 / 3))
$i = 3 * 5 Write-Host '3 * 5 = ' $i Write-Host '3 * 5 = ', $i Write-Host "3 * 5 = $i" Write-Output "3 * 5 = $i" Echo "3 * 5 = $i"
foreach ($i in 1..9) { Write-Host "$i, " -nonewline }
foreach ($i in 1..9) { if ($i % 3 -eq 0) { Write-Host "$i, " -nonewline } }
$sum = 0 foreach ($i in 1..99) { if ($i % 3 -eq 0) { $sum += $i } } Write-Host $sum $sum = 0 foreach ($i in 1..99) { $sum += if ($i % 3 -eq 0) {$i} else {0} } Write-Host $sum
1..9
1..9 | where { $_ % 3 -eq 0}
$sum = 0 1..99 | where {$_ % 3 -eq 0} | foreach {$sum += $_} $sum
# 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 function sn($a, $lim) { $n = [Math]::Floor($lim / $a) # 項数:n = 上限:lim / 公差:a $l = $n * $a # 末項:l = 項数:n * 公差:a ($a + $l) * $n / 2 # 総和:sn = (初項:a + 末項:l) * 項数:n / 2 } # 3 の倍数の合計 Write-Host (sn 3 999)
# 10000 までの 自然数の和 # 項目数 n = 10000 $n = 10000 Write-Host($n * ($n + 1) / 2)
# 10000 までの 偶数の和 # 項目数 n = 5000 $n = 10000 / 2 Write-Host($n * ($n + 1))
# 10000 までの 奇数の和 # 項目数 n = 5000 $n = 10000 / 2 Write-Host([Math]::Pow($n, 2))
# 1000 までの 自然数の2乗の和 $n = 1000 Write-Host($n * ($n + 1) * (2 * $n + 1) / 6)
# 100 までの 自然数の3乗の和 $n = 100 Write-Host([Math]::Pow($n, 2) * [Math]::Pow(($n + 1), 2) / 4)
# 初項 2, 公比 3, 項数 10 の等比数列の和 $n = 10; $a = 2; $r = 3; Write-Host(($a * ([Math]::Pow($r, $n) - 1)) / ($r - 1))
$a = 5 # 初項 5 $d = 3 # 公差 3 $n = 10 # 項数 10 $p = 1 # 積 foreach ($i in 1..$n) { $m = $a + ($d * ($i - 1)) $p *= $m; } Write-Host $p
# 初項 5, 公差 3, 項数 10 の数列の積を表示する Write-Host (product 5 3 10) function product($m, $d, $n) { if ($n -eq 0) { 1 } else { $m * (product ($m + $d) $d ($n - 1)) } }
# 階乗を求める関数 function Fact($n) { if ($n -le 1) { 1 } else { $n * (Fact ($n - 1)) } } # 10の階乗 Write-Host (Fact(10)) Write-Host (10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1)
# 下降階乗冪 function FallingFact($x, $n) { if ($n -le 1) { $x } else { $x * (FallingFact ($x - 1) ($n - 1)) } } # 10 から 6 までの 総乗 Write-Host (FallingFact 10 5) Write-Host (10 * 9 * 8 * 7 * 6)
# 上昇階乗冪 function RisingFact($x, $n) { if ($n -le 1) { $x } else { $x * (RisingFact ($x + 1) ($n - 1)) } } # 10 から 14 までの 総乗 Write-Host (RisingFact 10 5) Write-Host (10 * 11 * 12 * 13 * 14)
# 階乗 function Fact($n) { if ($n -le 1) { 1 } else { $n * (Fact ($n - 1)) } } # 下降階乗冪 function FallingFact($x, $n) { if ($n -le 1) { $x } else { $x * (FallingFact ($x - 1) ($n - 1)) } } # 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) $n = 10 $r = 5 Write-Host ((Fact $n) / (Fact ($n - $r))) Write-Host (FallingFact $n $r)
# 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) $n = 10 $r = 5 Write-Host ([Math]::Pow($n, $r))
# 組合せ function Comb($n, $r) { if (($r -eq 0) -or ($r -eq $n)) { 1 } elseif ($r -eq 1) { $n } else { (Comb ($n - 1) ($r - 1)) + (Comb ($n - 1) $r) } } # 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数) $n = 10 $r = 5 Write-Host (Comb $n $r)
# 組合せ function Comb($n, $r) { if (($r -eq 0) -or ($r -eq $n)) { 1 } elseif ($r -eq 1) { $n } else { (Comb ($n - 1) ($r - 1)) + (Comb ($n - 1) $r) } } # 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数) $n = 10 $r = 5 Write-Host (Comb ($n + $r - 1) $r)
# 自作の正弦関数 function mySin($x, $n, $nega, $numerator, $denominator, $y) { $m = 2 * $n $denominator = $denominator * ($m + 1) * $m $numerator = $numerator * $x * $x $a = $numerator / $denominator # 十分な精度になったら処理を抜ける if ($a -le 0.00000000001) { $y } else { if (-not $nega) { $a = -$a} $y + (mySin $x ($n + 1) (-not $nega) $numerator $denominator $a) } } foreach ($i in 0..24) { $degree = $i * 15 if ($degree % 30 -eq 0 -or $degree % 45 -eq 0) { $radian = $degree * [Math]::PI / 180.0 # 自作の正弦関数 $d1 = mySin $radian 1 $false $radian 1.0 $radian # 標準の正弦関数 $d2 = [Math]::Sin($radian) # 標準関数との差異 Write-Host ([String]::Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", $degree, $d1, $d2, $d1 - $d2)) } }
# 自作の余弦関数 function myCos($x, $n, $nega, $numerator, $denominator, $y) { $m = 2 * $n $denominator = $denominator * $m * ($m - 1) $numerator = $numerator * $x * $x $a = $numerator / $denominator # 十分な精度になったら処理を抜ける if ($a -le 0.00000000001) { $y } else { if (-not $nega) { $a = -$a} $y + (myCos $x ($n + 1) (-not $nega) $numerator $denominator $a) } } foreach ($i in 0..24) { $degree = $i * 15 if ($degree % 30 -eq 0 -or $degree % 45 -eq 0) { $radian = $degree * [Math]::PI / 180.0 # 自作の余弦関数 $d1 = myCos $radian 1 $false 1.0 1.0 1.0 # 標準の余弦関数 $d2 = [Math]::Cos($radian) # 標準関数との差異 Write-Host ([String]::Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", $degree, $d1, $d2, $d1 - $d2)) } }
# 自作の正接関数 function myTan($x, $x2, $n, $t) { $t = $x2 / ($n - $t) $n -= 2 if ($n -le 1) { $x / (1 - $t) } else { myTan $x $x2 $n $t } } foreach ($i in 0..12) { if (($i * 15) % 180 -ne 0) { $degree = $i * 15 - 90 $radian = $degree * [Math]::PI / 180.0 $x2 = $radian * $radian # 自作の正接関数 $d1 = myTan $radian $x2 15 0.0 # 15:必要な精度が得られる十分大きな奇数 # 標準の正接関数 $d2 = [Math]::Tan($radian) # 標準関数との差異 Write-Host ([String]::Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", $degree, $d1, $d2, $d1 - $d2)) } }
# 自作の指数関数 function myExp($x, $n, $numerator, $denominator, $y) { $denominator = $denominator * $n $numerator = $numerator * $x $a = $numerator / $denominator # 十分な精度になったら処理を抜ける if ([Math]::Abs($a) -le 0.00000000001) { $y } else { $y + (myExp $x ($n + 1) $numerator $denominator $a) } } foreach ($i in 0..20) { $x = ($i - 10) / 4.0 # 標準の指数関数 $d1 = [Math]::Exp($x) # 自作の指数関数 $d2 = myExp $x 1 1.0 1.0 1.0 # 標準関数との差異 Write-Host ([String]::Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", $x, $d1, $d2, $d1 - $d2)) }
# 自作の指数関数 function myExp($x, $x2, $n, $t) { $t = $x2 / ($n + $t) $n -= 4 if ($n -lt 6) { 1 + ((2 * $x) / (2 - $x + $t)) } else { myExp $x $x2 $n $t } } foreach ($i in 0..20) { $x = ($i - 10) / 4.0 # 標準の指数関数 $d1 = [Math]::Exp($x) # 自作の指数関数 $x2 = $x * $x $d2 = myExp $x $x2 30 0.0 # 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる # 標準関数との差異 Write-Host ([String]::Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", $x, $d1, $d2, $d1 - $d2)) }
# 自作の対数関数 function myLog($x2, $numerator, $denominator, $y) { $denominator = $denominator + 2 $numerator = $numerator * $x2 * $x2 $a = $numerator / $denominator # 十分な精度になったら処理を抜ける if ([Math]::Abs($a) -le 0.00000000001) { $y } else { $y + (myLog $x2 $numerator $denominator $a) } } foreach ($i in 1..20) { $x = $i / 5.0 # 標準の対数関数 $d1 = [Math]::Log($x) # 自作の対数関数 $x2 = ($x - 1) / ($x + 1) $d2 = 2 * (myLog $x2 $x2 1.0 $x2) # 標準関数との差異 Write-Host ([String]::Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", $x, $d1, $d2, $d1 - $d2)) }
# 自作の対数関数 function myLog($x, $n, $t) { $n2 = $n $x2 = $x if ($n -gt 3) { if ($n % 2 -eq 0) { $n2 = 2 } $x2 = $x * [Math]::Floor($n / 2) } $t = $x2 / ($n2 + $t) if ($n -le 2) { $x / (1 + $t) } else { myLog $x ($n - 1) $t } } foreach ($i in 1..20) { $x = $i / 5.0 # 標準の対数関数 $d1 = [Math]::Log($x) # 自作の対数関数 $d2 = myLog ($x - 1) 27 0.0 # 27:必要な精度が得られる十分大きな奇数 # 標準関数との差異 Write-Host ([String]::Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", $x, $d1, $d2, $d1 - $d2)) }
# 自作の双曲線正弦関数 function mySinh($x, $n, $numerator, $denominator, $y) { $m = 2 * $n $denominator = $denominator * ($m + 1) * $m $numerator = $numerator * $x * $x $a = $numerator / $denominator # 十分な精度になったら処理を抜ける if ([Math]::Abs($a) -le 0.00000000001) { $y } else { $y + (mySinh $x ($n + 1) $numerator $denominator $a) } } foreach ($i in 0..20) { $x = $i - 10 # 自作の双曲線正弦関数 $d1 = mySinh $x 1 $x 1.0 $x # 標準の双曲線正弦関数 $d2 = [Math]::Sinh($x) # 標準関数との差異 Write-Host ([String]::Format("{0,3:D} : {1,17:F10} - {2,17:F10} = {3,13:F10}", $x, $d1, $d2, $d1 - $d2)) }
# 自作の双曲線余弦関数 function myCosh($x, $n, $numerator, $denominator, $y) { $m = 2 * $n $denominator = $denominator * $m * ($m - 1) $numerator = $numerator * $x * $x $a = $numerator / $denominator # 十分な精度になったら処理を抜ける if ([Math]::Abs($a) -le 0.00000000001) { $y } else { $y + (myCosh $x ($n + 1) $numerator $denominator $a) } } foreach ($i in 0..20) { $x = $i - 10 # 自作の双曲線余弦関数 $d1 = myCosh $x 1 1.0 1.0 1.0 # 標準の双曲線余弦関数 $d2 = [Math]::Cosh($x) # 標準関数との差異 Write-Host ([String]::Format("{0,3:D} : {1,17:F10} - {2,17:F10} = {3,13:F10}", $x, $d1, $d2, $d1 - $d2)) }
function f($x) { 4 / (1 + $x * $x) } $a = 0 $b = 1 # 台形則で積分 $n = 2; foreach ($j in 1..10) { $h = ($b - $a) / $n $s = 0 $x = $a foreach ($i in 1..($n - 1)) { $x += $h $s += (f $x) } $s = $h * (((f $a) + (f $b)) / 2 + $s) $n *= 2 # 結果を π と比較 Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, $s - [Math]::PI)) }
function f($x) { 4 / (1 + $x * $x) } $a = 0 $b = 1 # 中点則で積分 $n = 2; foreach ($j in 1..10) { $h = ($b - $a) / $n $s = 0 $x = $a + ($h / 2) foreach ($i in 1..$n) { $s += (f $x) $x += $h } $s *= $h $n *= 2 # 結果を π と比較 Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, $s - [Math]::PI)) }
function f($x) { 4 / (1 + $x * $x) } $a = 0 $b = 1 # Simpson則で積分 $n = 2; foreach ($j in 1..5) { $h = ($b - $a) / $n $s2 = 0 $s4 = 0 $x = $a + $h foreach ($i in 1..($n / 2)) { $s4 += (f $x) $x += $h $s2 += (f $x) $x += $h } $s2 = ($s2 - (f $b)) * 2 + (f $a) + (f $b) $s4 *= 4 $s = ($s2 + $s4) * $h / 3 $n *= 2 # 結果を π と比較 Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, $s - [Math]::PI)) }
function f($x) { 4 / (1 + $x * $x) } $a = 0 $b = 1 $t = New-Object "double[,]" 7,7 # 台形則で積分 $n = 2; foreach ($i in 1..6) { $h = ($b - $a) / $n $s = 0 $x = $a foreach ($j in 1..($n - 1)) { $x += $h $s += (f $x) } # 結果を保存 $t[$i,1] = $h * (((f $a) + (f $b)) / 2 + $s) $n *= 2 } # Richardsonの補外法 $n = 4 foreach ($j in 2..6) { foreach ($i in $j..6) { $t1 = $t[ $i ,($j - 1)] $t2 = $t[($i - 1),($j - 1)] $t[$i,$j] = $t1 + ($t1 - $t2) / ($n - 1) if ($i -eq $j) { # 結果を π と比較 $s = $t[$i,$j] Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, ($s - [Math]::PI))) } } $n *= 4; }
# データ点の数 set-variable -option constant -name N -value 7 # 元の関数 function f($x) { $x - [Math]::Pow($x, 3) / (3 * 2) + [Math]::Pow($x, 5) / (5 * 4 * 3 * 2) } # Lagrange (ラグランジュ) 補間 function lagrange($d, $x, $y) { $sum = 0 foreach ($i in 0..($N - 1)) { $prod = $y[$i] foreach ($j in 0..($N - 1)) { if ($j -ne $i) { $prod *= ($d - $x[$j]) / ($x[$i] - $x[$j]) } } $sum += $prod } $sum } $x = New-Object double[] $N $y = New-Object double[] $N # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット foreach ($i in 0..($N - 1)) { $d = $i * 1.5 - 4.5 $x[$i] = $d $y[$i] = f($d) } # 0.5刻みで 与えられていない値を補間 foreach ($i in 0..18) { $d = $i * 0.5 - 4.5 $d1 = f($d) $d2 = (lagrange $d $x $y) # 元の関数と比較 Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d, $d1, $d2, ($d1 - $d2))) }
# データ点の数 set-variable -option constant -name N -value 7 # 元の関数 function f($x) { $x - [Math]::Pow($x, 3) / (3 * 2) + [Math]::Pow($x, 5) / (5 * 4 * 3 * 2) } # Neville (ネヴィル) 補間 function neville($d, $x, $y) { $w = New-Object "double[,]" $N,$N foreach ($i in 0..($N - 1)) { $w[0,$i] = $y[$i] } foreach ($j in 1..($N - 1)) { foreach ($i in 0..($N - $j - 1)) { $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-1),0] } $x = New-Object double[] $N $y = New-Object double[] $N # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット foreach ($i in 0..($N - 1)) { $d = $i * 1.5 - 4.5 $x[$i] = $d $y[$i] = f($d) } # 0.5刻みで 与えられていない値を補間 foreach ($i in 0..18) { $d = $i * 0.5 - 4.5 $d1 = f($d) $d2 = (neville $d $x $y) # 元の関数と比較 Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d, $d1, $d2, ($d1 - $d2))) }
# データ点の数 set-variable -option constant -name N -value 7 # 元の関数 function f($x) { $x - [Math]::Pow($x, 3) / (3 * 2) + [Math]::Pow($x, 5) / (5 * 4 * 3 * 2) } # Newton (ニュートン) 補間 function newton($d, $x, $a) { $sum = $a[0] foreach ($i in 1..($N - 1)) { $prod = $a[$i] foreach ($j in 0..($i - 1)) { $prod *= ($d - $x[$j]) } $sum += $prod } $sum } $x = New-Object double[] $N $y = New-Object double[] $N # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット foreach ($i in 0..($N - 1)) { $d1 = $i * 1.5 - 4.5 $x[$i] = $d1 $y[$i] = f($d1) } # 差分商の表を作る $d = New-Object "double[,]" $N,$N foreach ($j in 0..($N - 1)) { $d[0,$j] = $y[$j] } foreach ($i in 1..($N - 1)) { foreach ($j in 0..($N - $i - 1)) { $d[$i,$j] = ($d[($i-1),($j+1)] - $d[($i-1),$j]) / ($x[($j+$i)] - $x[$j]) } } # n階差分商 $a = New-Object double[] $N foreach ($j in 0..($N - 1)) { $a[$j] = $d[$j,0] } # 0.5刻みで 与えられていない値を補間 foreach ($i in 0..18) { $d1 = $i * 0.5 - 4.5 $d2 = f($d1) $d3 = (newton $d1 $x $a) # 元の関数と比較 Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d1, $d2, $d3, ($d2 - $d3))) }
# データ点の数 set-variable -option constant -name N -value 7 set-variable -option constant -name Nx2 -value 14 # 元の関数 function f($x) { $x - [Math]::Pow($x, 3) / (3 * 2) + [Math]::Pow($x, 5) / (5 * 4 * 3 * 2) } # 導関数 function fd($x) { 1 - [Math]::Pow($x, 2) / 2 + [Math]::Pow($x, 4) / (4 * 3 * 2) } # Hermite (エルミート) 補間 function hermite($d, $z, $a) { $sum = $a[0] foreach ($i in 1..($Nx2 - 1)) { $prod = $a[$i] foreach ($j in 0..($i - 1)) { $prod *= ($d - $z[$j]) } $sum += $prod } $sum } $x = New-Object double[] $N $y = New-Object double[] $N $yd = New-Object double[] $N # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット foreach ($i in 0..($N - 1)) { $d1 = $i * 1.5 - 4.5 $x[$i] = $d1 $y[$i] = f($d1) $yd[$i] = fd($d1) } # 差分商の表を作る $z = New-Object double[] $Nx2 $d = New-Object "double[,]" $Nx2,$Nx2 foreach ($i in 0..($Nx2 - 1)) { $j = [Math]::Floor($i / 2) $z[$i] = $x[$j] $d[0,$i] = $y[$j] } foreach ($i in 1..($Nx2 - 1)) { foreach ($j in 0..($Nx2 - $i - 1)) { if ($i -eq 1 -and $j % 2 -eq 0) { $d[$i,$j] = $yd[[Math]::Floor($j / 2)] } else { $d[$i,$j] = ($d[($i-1),($j+1)] - $d[($i-1),$j]) / ($z[($j+$i)] - $z[$j]) } } } # n階差分商 $a = New-Object double[] $Nx2 foreach ($j in 0..($Nx2 - 1)) { $a[$j] = $d[$j,0] } # 0.5刻みで 与えられていない値を補間 foreach ($i in 0..18) { $d1 = $i * 0.5 - 4.5 $d2 = f($d1) $d3 = (hermite $d1 $z $a) # 元の関数と比較 Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d1, $d2, $d3, ($d2 - $d3))) }
# データ点の数 set-variable -option constant -name N -value 7 # 元の関数 function f($x) { $x - [Math]::Pow($x, 3) / (3 * 2) + [Math]::Pow($x, 5) / (5 * 4 * 3 * 2) } # Spline (スプライン) 補間 function spline($d, $x, $y, $z) { # 補間関数値がどの区間にあるか $k = -1 foreach ($i in 1..($N - 1)) { if ($d -le $x[$i]) { $k = $i - 1 break } } if ($k -lt 0) { $k = $N - 1 } $d1 = $x[($k+1)] - $d $d2 = $d - $x[$k] $d3 = $x[($k+1)] - $x[$k] ($z[$k] * [Math]::Pow($d1,3) + $z[($k+1)] * [Math]::Pow($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 } $x = New-Object double[] $N $y = New-Object double[] $N # 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット foreach ($i in 0..($N - 1)) { $d1 = $i * 1.5 - 4.5 $x[$i] = $d1 $y[$i] = f($d1) } # 3項方程式の係数の表を作る $a = New-Object double[] $N $b = New-Object double[] $N $c = New-Object double[] $N $d = New-Object double[] $N foreach ($i in 1..($N - 2)) { $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項方程式を解く (ト−マス法) $g = New-Object double[] $N $s = New-Object double[] $N $g[1] = $b[1] $s[1] = $d[1] foreach ($i in 2..($N - 2)) { $g[$i] = $b[$i] - $a[$i] * $c[($i-1)] / $g[($i-1)] $s[$i] = $d[$i] - $a[$i] * $s[($i-1)] / $g[($i-1)] } $z = New-Object double[] $N $z[0] = 0 $z[($N-1)] = 0 $z[($N-2)] = $s[($N-2)] / $g[($N-2)] for ($i = $N - 3; $i -ge 1; $i--) { $z[$i] = ($s[$i] - $c[$i] * $z[($i+1)]) / $g[$i] } # 0.5刻みで 与えられていない値を補間 foreach ($i in 0..18) { $d1 = $i * 0.5 - 4.5 $d2 = f($d1) $d3 = (spline $d1 $x $y $z) # 元の関数と比較 Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d1, $d2, $d3, ($d2 - $d3))) }
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) $h = 0.01 # 角度 $degree = 45 $radian = $degree * [Math]::PI / 180.0 # 初速 250 km/h -> 秒速に変換 $v = [Math]::Floor(250 * 1000 / 3600) # 水平方向の速度 $vx = New-Object double[] 2 $vx[0] = $v * [Math]::Cos($radian) # 鉛直方向の速度 $vy = New-Object double[] 2 $vy[0] = $v * [Math]::Sin($radian) # 経過秒数 $t = 0.0 # 位置 $x = 0.0 $y = 0.0 # 空気抵抗による水平方向の減速分 function fx($vx, $vy) { return $global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vx } # 重力と空気抵抗による鉛直方向の減速分 function fy($vx, $vy) { return $global:g + ($global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vy) } # Euler法 for ($i = 1; $y -ge 0.0; $i++) { # 経過秒数 $t = $i * $h # 位置 $x += $h * $vx[0] $y += $h * $vy[0] Write-Host ([String]::Format("{0,4:F2}`t{1,8:F5}`t{2,9:F5}`t{3,9:F5}`t{4,8:F5}", $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] }
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) $h = 0.01 # 角度 $degree = 45 $radian = $degree * [Math]::PI / 180.0 # 初速 250 km/h -> 秒速に変換 $v = [Math]::Floor(250 * 1000 / 3600) # 水平方向の速度 $vx = New-Object double[] 3 $vx[0] = $v * [Math]::Cos($radian) # 鉛直方向の速度 $vy = New-Object double[] 3 $vy[0] = $v * [Math]::Sin($radian) # 経過秒数 $t = 0.0 # 位置 $x = New-Object double[] 3 $y = New-Object double[] 3 $x[0] = 0.0 $y[0] = 0.0 # 空気抵抗による水平方向の減速分 function fx($vx, $vy) { return $global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vx } # 重力と空気抵抗による鉛直方向の減速分 function fy($vx, $vy) { return $global:g + ($global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vy) } # Heun法 for ($i = 1; $y[0] -ge 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] Write-Host ([String]::Format("{0,4:F2}`t{1,8:F5}`t{2,9:F5}`t{3,9:F5}`t{4,8:F5}", $t, $vx[0], $vy[0], $x[0], $y[0])) }
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) $h = 0.01 # 角度 $degree = 45 $radian = $degree * [Math]::PI / 180.0 # 初速 250 km/h -> 秒速に変換 $v = [Math]::Floor(250 * 1000 / 3600) # 水平方向の速度 $vx = New-Object double[] 2 $vx[0] = $v * [Math]::Cos($radian) # 鉛直方向の速度 $vy = New-Object double[] 2 $vy[0] = $v * [Math]::Sin($radian) # 経過秒数 $t = 0.0 # 位置 $x = New-Object double[] 2 $y = New-Object double[] 2 $x[0] = 0.0 $y[0] = 0.0 # 空気抵抗による水平方向の減速分 function fx($vx, $vy) { return $global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vx } # 重力と空気抵抗による鉛直方向の減速分 function fy($vx, $vy) { return $global:g + ($global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vy) } # 中点法 for ($i = 1; $y[0] -ge 0.0; $i++) { # 経過秒数 $t = $i * $h # 位置・速度 $vx[1] = $h * (fx $vx[0] $vy[0]) $vy[1] = $h * (fy $vx[0] $vy[0]) $wx = $vx[0] + $vx[1] / 2 $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 Write-Host ([String]::Format("{0,4:F2}`t{1,8:F5}`t{2,9:F5}`t{3,9:F5}`t{4,8:F5}", $t, $vx[0], $vy[0], $x[0], $y[0])) }
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) $h = 0.01 # 角度 $degree = 45 $radian = $degree * [Math]::PI / 180.0 # 初速 250 km/h -> 秒速に変換 $v = [Math]::Floor(250 * 1000 / 3600) # 水平方向の速度 $vx = New-Object double[] 5 $vx[0] = $v * [Math]::Cos($radian) # 鉛直方向の速度 $vy = New-Object double[] 5 $vy[0] = $v * [Math]::Sin($radian) # 経過秒数 $t = 0.0 # 位置 $x = New-Object double[] 5 $y = New-Object double[] 5 $x[0] = 0.0 $y[0] = 0.0 # 空気抵抗による水平方向の減速分 function fx($vx, $vy) { return $global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vx } # 重力と空気抵抗による鉛直方向の減速分 function fy($vx, $vy) { return $global:g + ($global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vy) } # Runge-Kutta法 for ($i = 1; $y[0] -ge 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]) $wx = $vx[0] + $vx[1] / 2 $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 Write-Host ([String]::Format("{0,4:F2}`t{1,8:F5}`t{2,9:F5}`t{3,9:F5}`t{4,8:F5}", $t, $vx[0], $vy[0], $x[0], $y[0])) }
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) $h = 0.01 # 角度 $degree = 45 $radian = $degree * [Math]::PI / 180.0 # 初速 250 km/h -> 秒速に変換 $v = [Math]::Floor(250 * 1000 / 3600) # 水平方向の速度 $vx = New-Object double[] 5 $vx[0] = $v * [Math]::Cos($radian) # 鉛直方向の速度 $vy = New-Object double[] 5 $vy[0] = $v * [Math]::Sin($radian) # 経過秒数 $t = 0.0 # 位置 $x = New-Object double[] 5 $y = New-Object double[] 5 $x[0] = 0.0 $y[0] = 0.0 # 空気抵抗による水平方向の減速分 function fx($vx, $vy) { return $global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vx } # 重力と空気抵抗による鉛直方向の減速分 function fy($vx, $vy) { return $global:g + ($global:k * [Math]::Sqrt($vx * $vx + $vy * $vy) * $vy) } # Runge-Kutta-Gill法 for ($i = 1; $y[0] -ge 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]) $wx = $vx[0] + $vx[1] / 2 $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] * (([Math]::Sqrt(2.0) - 1) / 2) + $vx[2] * (1 - 1 / [Math]::Sqrt(2.0)) $wy = $vy[0] + $vy[1] * (([Math]::Sqrt(2.0) - 1) / 2) + $vy[2] * (1 - 1 / [Math]::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] / [Math]::Sqrt(2.0) + $vx[3] * (1 + 1 / [Math]::Sqrt(2.0)) $wy = $vy[0] - $vy[2] / [Math]::Sqrt(2.0) + $vy[3] * (1 + 1 / [Math]::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 - [Math]::Sqrt(2.0)) + $x[3] * (2 + [Math]::Sqrt(2.0)) + $x[4]) / 6 $y[0] += ( $y[1] + $y[2] * (2 - [Math]::Sqrt(2.0)) + $y[3] * (2 + [Math]::Sqrt(2.0)) + $y[4]) / 6 $vx[0] += ($vx[1] + $vx[2] * (2 - [Math]::Sqrt(2.0)) + $vx[3] * (2 + [Math]::Sqrt(2.0)) + $vx[4]) / 6 $vy[0] += ($vy[1] + $vy[2] * (2 - [Math]::Sqrt(2.0)) + $vy[3] * (2 + [Math]::Sqrt(2.0)) + $vy[4]) / 6 Write-Host ([String]::Format("{0,4:F2}`t{1,8:F5}`t{2,9:F5}`t{3,9:F5}`t{4,8:F5}", $t, $vx[0], $vy[0], $x[0], $y[0])) }
function bisection($a, $b) { while ($true) { # 区間 (a, b) の中点 c = (a + b) / 2 $c = ($a + $b) / 2 Write-Host ([String]::Format("{0,12:F10}`t{1,13:F10}", $c, $c - [Math]::Sqrt(2))) $fc = f($c) if ([Math]::Abs($fc) -lt 0.0000000001) { break } if ($fc -lt 0) { # f(c) < 0 であれば, 解は区間 (c, b) の中に存在 $a = $c } else { # f(c) > 0 であれば, 解は区間 (a, c) の中に存在 $b = $c } } return $c } function f($x) { return $x * $x - 2 } $a = 1 $b = 2 Write-Host ([String]::Format("{0,12:F10}", (bisection $a $b)))
function falseposition($a, $b) { while ($true) { # 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 $c = ($a * (f $b) - $b * (f $a)) / ((f $b) - (f $a)) Write-Host ([String]::Format("{0,12:F10}`t{1,13:F10}", $c, $c - [Math]::Sqrt(2))) $fc = f($c) if ([Math]::Abs($fc) -lt 0.0000000001) { break } if ($fc -lt 0) { # f(c) < 0 であれば, 解は区間 (c, b) の中に存在 $a = $c } else { # f(c) > 0 であれば, 解は区間 (a, c) の中に存在 $b = $c } } return $c } function f($x) { return $x * $x - 2 } $a = 1 $b = 2 Write-Host ([String]::Format("{0,12:F10}", (falseposition $a $b)))
function iterative($x0) { while ($true) { $x1 = g($x0) Write-Host ([String]::Format("{0,12:F10}`t{1,13:F10}", $x1, $x1 - [Math]::Sqrt(2))) if ([Math]::Abs($x1 - $x0) -lt 0.0000000001) { break } $x0 = $x1 } return $x1 } function g($x) { return ($x / 2) + (1 / $x) } $x = 1 Write-Host ([String]::Format("{0,12:F10}", (iterative $x)))
function newton($x0) { while ($true) { $x1 = $x0 - ((f0 $x0) / (f1 $x0)) Write-Host ([String]::Format("{0,12:F10}`t{1,13:F10}", $x1, $x1 - [Math]::Sqrt(2))) if ([Math]::Abs($x1 - $x0) -lt 0.0000000001) { break } $x0 = $x1 } return $x1 } function f0($x) { return $x * $x - 2 } function f1($x) { return 2 * $x } $x = 2 Write-Host ([String]::Format("{0,12:F10}", (newton $x)))
function bailey($x0) { while ($true) { $x1 = $x0 - ((f0 $x0) / ((f1 $x0) - ((f0 $x0) * (f2 $x0) / (2 * (f1 $x0))))) Write-Host ([String]::Format("{0,12:F10}`t{1,13:F10}", $x1, $x1 - [Math]::Sqrt(2))) if ([Math]::Abs($x1 - $x0) -lt 0.0000000001) { break } $x0 = $x1 } return $x1 } function f0($x) { return $x * $x - 2 } function f1($x) { return 2 * $x } function f2($x) { return 2 } $x = 2 Write-Host ([String]::Format("{0,12:F10}", (bailey $x)))
function secant($x0, $x1) { while ($true) { $x2 = $x1 - (f $x1) * ($x1 - $x0) / ((f $x1) - (f $x0)) Write-Host ([String]::Format("{0,12:F10}`t{1,13:F10}", $x2, $x2 - [Math]::Sqrt(2))) if ([Math]::Abs($x2 - $x1) -lt 0.0000000001) { break } $x0 = $x1 $x1 = $x2 } return $x2 } function f($x) { return $x * $x - 2 } $x0 = 1 $x1 = 2 Write-Host ([String]::Format("{0,12:F10}", (secant $x0 $x1)))
set-variable -name N -value 3 -option constant # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # ヤコビの反復法 function jacobi($a, $b, $x0) { while ($true) { $x1 = 0.0, 0.0, 0.0, 0.0 $finish = $true foreach ($i in 0..$N) { $x1[$i] = 0 foreach ($j in 0..$N) { if ($j -ne $i) { $x1[$i] += $a[$i][$j] * $x0[$j] } } $x1[$i] = ($b[$i] - $x1[$i]) / $a[$i][$i] if ([Math]::Abs($x1[$i] - $x0[$i]) -gt 0.0000000001) { $finish = $false } } foreach ($i in 0..$N) { $x0[$i] = $x1[$i] } if ($finish) { return } disp_vector $x0 } } $a = (9.0,2.0,1.0,1.0), (2.0,8.0,-2.0,1.0), (-1.0,-2.0,7.0,-2.0), (1.0,-1.0,-2.0,6.0) $b = 20.0, 16.0, 8.0, 17.0 $c = 0.0, 0.0, 0.0, 0.0 # ヤコビの反復法 jacobi $a $b $c Write-Host "X" disp_vector $c
set-variable -name N -value 3 -option constant # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # ガウス・ザイデル法 function gauss($a, $b, $x0) { while ($true) { $finish = $true foreach ($i in 0..$N) { $x1 = 0 foreach ($j in 0..$N) { if ($j -ne $i) { $x1 += $a[$i][$j] * $x0[$j] } } $x1 = ($b[$i] - $x1) / $a[$i][$i] if ([Math]::Abs($x1 - $x0[$i]) -gt 0.0000000001) { $finish = $false } $x0[$i] = $x1 } if ($finish) { return } disp_vector $x0 } } $a = (9.0,2.0,1.0,1.0), (2.0,8.0,-2.0,1.0), (-1.0,-2.0,7.0,-2.0), (1.0,-1.0,-2.0,6.0) $b = 20.0, 16.0, 8.0, 17.0 $c = 0.0, 0.0, 0.0, 0.0 # ガウス・ザイデル法 gauss $a $b $c Write-Host "X" disp_vector $c
set-variable -name N -value 3 -option constant # 前進消去 function forward_elimination($a, $b) { foreach ($pivot in 0..($N - 1)) { foreach ($row in ($pivot + 1)..$N) { $s = $a[$row][$pivot] / $a[$pivot][$pivot] foreach ($col in $pivot..$N) { $a[$row][$col] -= $a[$pivot][$col] * $s } $b[$row] -= $b[$pivot] * $s } } } # 後退代入 function backward_substitution($a, $b) { foreach ($row in $N..0) { if (($row + 1) -le $N) { foreach ($col in $N..($row + 1)) { $b[$row] -= $a[$row][$col] * $b[$col] } } $b[$row] /= $a[$row][$row] } } # ピボット選択 function pivoting($a, $b) { foreach ($pivot in 0..$N) { # 各列で 一番値が大きい行を 探す $max_row = $pivot $max_val = 0 foreach ($row in $pivot..$N) { if ([Math]::Abs($a[$row][$pivot]) -gt $max_val) { # 一番値が大きい行 $max_val = [Math]::Abs($a[$row][$pivot]) $max_row = $row } } # 一番値が大きい行と入れ替え if ($max_row -ne $pivot) { foreach ($col in 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 } } } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } $a = (-1.0,-2.0,7.0,-2.0), (1.0,-1.0,-2.0,6.0), (9.0,2.0,1.0,1.0), (2.0,8.0,-2.0,1.0) $b = 8.0, 17.0, 20.0, 16.0 # ピボット選択 pivoting $a $b Write-Host "pivoting" Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 前進消去 forward_elimination $a $b Write-Host "forward elimination" Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 後退代入 backward_substitution $a $b Write-Host "X" disp_vector $b
set-variable -name N -value 3 -option constant # 前進消去 function forward_elimination($a, $b) { foreach ($pivot in 0..$N) { foreach ($row in 0..$N) { if ($row -eq $pivot) { continue } $s = $a[$row][$pivot] / $a[$pivot][$pivot] foreach ($col in $pivot..$N) { $a[$row][$col] -= $a[$pivot][$col] * $s } $b[$row] -= $b[$pivot] * $s } } } # 後退代入 function backward_substitution($a, $b) { foreach ($pivot in 0..$N) { $b[$pivot] /= $a[$pivot][$pivot] } } # ピボット選択 function pivoting($a, $b) { foreach ($pivot in 0..$N) { # 各列で 一番値が大きい行を 探す $max_row = $pivot $max_val = 0 foreach ($row in $pivot..$N) { if ([Math]::Abs($a[$row][$pivot]) -gt $max_val) { # 一番値が大きい行 $max_val = [Math]::Abs($a[$row][$pivot]) $max_row = $row } } # 一番値が大きい行と入れ替え if ($max_row -ne $pivot) { foreach ($col in 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 } } } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } $a = (-1.0,-2.0,7.0,-2.0), (1.0,-1.0,-2.0,6.0), (9.0,2.0,1.0,1.0), (2.0,8.0,-2.0,1.0) $b = 8.0, 17.0, 20.0, 16.0 # ピボット選択 pivoting $a $b Write-Host "pivoting" Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 前進消去 forward_elimination $a $b Write-Host "forward elimination" Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 後退代入 backward_substitution $a $b Write-Host "X" disp_vector $b
set-variable -name N -value 3 -option constant # 前進消去 function forward_elimination($a, $b) { foreach ($pivot in 0..($N - 1)) { foreach ($row in ($pivot + 1)..$N) { $s = $a[$row][$pivot] / $a[$pivot][$pivot] foreach ($col in $pivot..$N) { $a[$row][$col] -= $a[$pivot][$col] * $s # これが 上三角行列 } $a[$row][$pivot] = $s # これが 下三角行列 # $b[$row] -= $b[$pivot] * $s # この値は変更しない } } } # 前進代入 function forward_substitution($a, $b, $y) { foreach ($row in 0..$N) { foreach ($col in 0..$row) { $b[$row] -= $a[$row][$col] * $y[$col] } $y[$row] = $b[$row] } } # 後退代入 function backward_substitution($a, $y, $x) { foreach ($row in $N..0) { if (($row+1) -le $N) { foreach ($col in $N..($row + 1)) { $y[$row] -= $a[$row][$col] * $x[$col] } } $x[$row] = $y[$row] / $a[$row][$row] } } # ピボット選択 function pivoting($a, $b) { foreach ($pivot in 0..$N) { # 各列で 一番値が大きい行を 探す $max_row = $pivot $max_val = 0 foreach ($row in $pivot..$N) { if ([Math]::Abs($a[$row][$pivot]) -gt $max_val) { # 一番値が大きい行 $max_val = [Math]::Abs($a[$row][$pivot]) $max_row = $row } } # 一番値が大きい行と入れ替え if ($max_row -ne $pivot) { foreach ($col in 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 } } } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } $a = (-1.0,-2.0,7.0,-2.0), (1.0,-1.0,-2.0,6.0), (9.0,2.0,1.0,1.0), (2.0,8.0,-2.0,1.0) $b = 8.0, 17.0, 20.0, 16.0 # ピボット選択 pivoting $a $b Write-Host "pivoting" Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 前進消去 forward_elimination $a $b Write-Host "LU" disp_matrix $a # Ly=b から y を求める (前進代入) $y = 0.0, 0.0, 0.0, 0.0 forward_substitution $a $b $y Write-Host "Y" disp_vector $y # Ux=y から x を求める (後退代入) $x = 0.0, 0.0, 0.0, 0.0 backward_substitution $a $y $x Write-Host "X" disp_vector $x
set-variable -name N -value 3 -option constant # 前進消去 function forward_elimination($a, $b) { foreach ($pivot in 0..$N) { $s = 0 for ($col = 0; $col -lt $pivot; $col++) { $s += $a[$pivot][$col] * $a[$pivot][$col] } # ここで根号の中が負の値になると計算できない! $a[$pivot][$pivot] = [Math]::Sqrt($a[$pivot][$pivot] - $s) for ($row = $pivot + 1; $row -le $N; $row++) { $s = 0 for ($col = 0; $col -lt $pivot; $col++) { $s += $a[$row][$col] * $a[$pivot][$col] } $a[$row][$pivot] = ($a[$row][$pivot] - $s) / $a[$pivot][$pivot] $a[$pivot][$row] = $a[$row][$pivot] } } } # 前進代入 function forward_substitution($a, $b, $y) { foreach ($row in 0..$N) { foreach ($col in 0..($row - 1)) { $b[$row] -= $a[$row][$col] * $y[$col] } $y[$row] = $b[$row] / $a[$row][$row] } } # 後退代入 function backward_substitution($a, $y, $x) { foreach ($row in $N..0) { if (($row+1) -le $N) { foreach ($col in $N..($row + 1)) { $y[$row] -= $a[$row][$col] * $x[$col] } } $x[$row] = $y[$row] / $a[$row][$row] } } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } $a = (5.0, 2.0, 3.0, 4.0), (2.0, 10.0, 6.0, 7.0), (3.0, 6.0, 15.0, 9.0), (4.0, 7.0, 9.0, 20.0) $b = 34.0, 68.0, 96.0, 125.0 Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 前進消去 forward_elimination $a $b Write-Host "LL^T" disp_matrix $a # Ly=b から y を求める (前進代入) $y = 0.0, 0.0, 0.0, 0.0 forward_substitution $a $b $y Write-Host "Y" disp_vector $y # L^Tx=y から x を求める (後退代入) $x = 0.0, 0.0, 0.0, 0.0 backward_substitution $a $y $x Write-Host "X" disp_vector $x
set-variable -name N -value 3 -option constant # 前進消去 function forward_elimination($a, $b) { foreach ($pivot in 0..$N) { # pivot < k の場合 $s = 0 for ($col = 0; $col -lt $pivot; $col++) { $s = $a[$pivot][$col] for ($k = 0; $k -lt $col; $k++) { $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; $k -lt $pivot; $k++) { foreach ($k in 0..($pivot-1)) { $s -= $a[$pivot][$k] * $a[$pivot][$k] * $a[$k][$k] } } $a[$pivot][$pivot] = $s } } # 前進代入 function forward_substitution($a, $b, $y) { foreach ($row in 0..$N) { foreach ($col in 0..($row - 1)) { $b[$row] -= $a[$row][$col] * $y[$col] } $y[$row] = $b[$row] } } # 後退代入 function backward_substitution($a, $y, $x) { foreach ($row in $N..0) { if (($row+1) -le $N) { foreach ($col in $N..($row + 1)) { $y[$row] -= $a[$row][$col] * $a[$row][$row] * $x[$col] } } $x[$row] = $y[$row] / $a[$row][$row] } } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } $a = (5.0, 2.0, 3.0, 4.0), (2.0, 10.0, 6.0, 7.0), (3.0, 6.0, 15.0, 9.0), (4.0, 7.0, 9.0, 20.0) $b = 34.0, 68.0, 96.0, 125.0 Write-Host "A" disp_matrix $a Write-Host "B" disp_vector $b Write-Host "" # 前進消去 forward_elimination $a $b Write-Host "LDL^T" disp_matrix $a # Ly=b から y を求める (前進代入) $y = 0.0, 0.0, 0.0, 0.0 forward_substitution $a $b $y Write-Host "Y" disp_vector $y # DL^Tx=y から x を求める (後退代入) $x = 0.0, 0.0, 0.0, 0.0 backward_substitution $a $y $x Write-Host "X" disp_vector $x
set-variable -name N -value 3 -option constant # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 正規化 (ベクトルの長さを1にする) function normarize($x) { $s = 0.0 foreach ($i in 0..$N) { $s += $x[$i] * $x[$i] } $s = [Math]::Sqrt($s) foreach ($i in 0..$N) { $x[$i] /= $s } } # ベキ乗法 function power($a, $x0) { $lambda = 0.0 # 正規化 (ベクトル x0 の長さを1にする) normarize $x0 $e0 = 0.0 foreach ($i in 0..$N) { $e0 += $x0[$i] } foreach ($k in 1..200) { # 1次元配列を表示 Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_vector $x0 # 行列の積 x1 = A × x0 $x1 = (0.0, 0.0, 0.0, 0.0) foreach ($i in 0..$N) { foreach ($j in 0..$N) { $x1[$i] += $a[$i][$j] * $x0[$j] } } # 内積 $p0 = 0.0 $p1 = 0.0 foreach ($i in 0..$N) { $p0 += $x1[$i] * $x1[$i] $p1 += $x1[$i] * $x0[$i] } # 固有値 $lambda = $p0 / $p1 # 正規化 (ベクトル x1 の長さを1にする) normarize $x1 # 収束判定 $e1 = 0.0 foreach ($i in 0..$N) { $e1 += $x1[$i] } if ([Math]::Abs($e0 - $e1) -lt 0.00000000001) { break } foreach ($i in 0..$N) { $x0[$i] = $x1[$i] } $e0 = $e1 } $lambda } # ベキ乗法で最大固有値を求める $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) $x = (1.0 ,0.0 ,0.0 ,0.0) # ベキ乗法 $lambda = (power $a $x) Write-Host Write-Host "eigenvalue" Write-Host ([String]::Format("{0,14:F10}", $lambda)) Write-Host "eigenvector" disp_vector $x
set-variable -name N -value 3 -option constant # LU分解 function forward_elimination($a) { foreach ($pivot in 0..($N - 1)) { foreach ($row in ($pivot + 1)..$N) { $s = $a[$row][$pivot] / $a[$pivot][$pivot] foreach ($col in $pivot..$N) { $a[$row][$col] -= $a[$pivot][$col] * $s # これが 上三角行列 } $a[$row][$pivot] = $s # これが 下三角行列 } } } # 前進代入 function forward_substitution($a, $y, $b) { foreach ($row in 0..$N) { foreach ($col in 0..$row) { $b[$row] -= $a[$row][$col] * $y[$col] } $y[$row] = $b[$row] } } # 後退代入 function backward_substitution($a, $x, $y) { foreach ($row in $N..0) { if ($row -lt $N) { foreach ($col in $N..($row + 1)) { $y[$row] -= $a[$row][$col] * $x[$col] } } $x[$row] = $y[$row] / $a[$row][$row] } } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 正規化 (ベクトルの長さを1にする) function normarize($x) { $s = 0.0 foreach ($i in 0..$N) { $s += $x[$i] * $x[$i] } $s = [Math]::Sqrt($s) foreach ($i in 0..$N) { $x[$i] /= $s } } # 逆ベキ乗法 function inverse($a, $x0) { $lambda = 0.0 # 正規化 (ベクトル x0 の長さを1にする) normarize $x0 $e0 = 0.0 foreach ($i in 0..$N) { $e0 += $x0[$i] } foreach ($k in 1..200) { # 1次元配列を表示 Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_vector $x0 # Ly = b から y を求める (前進代入) $y = (0.0, 0.0, 0.0, 0.0) $b = (0.0, 0.0, 0.0, 0.0) foreach ($i in 0..$N) { $b[$i] = $x0[$i] } forward_substitution $a $y $b # Ux = y から x を求める (後退代入) $x1 = (0.0, 0.0, 0.0, 0.0) backward_substitution $a $x1 $y # 内積 $p0 = 0.0 $p1 = 0.0 foreach ($i in 0..$N) { $p0 += $x1[$i] * $x1[$i] $p1 += $x1[$i] * $x0[$i] } # 固有値 $lambda = $p1 / $p0 # 正規化 (ベクトル x1 の長さを1にする) normarize $x1 # 収束判定 $e1 = 0.0 foreach ($i in 0..$N) { $e1 += $x1[$i] } if ([Math]::Abs($e0 - $e1) -lt 0.00000000001) { break } $lambda0 = $lambda1 foreach ($i in 0..$N) { $x0[$i] = $x1[$i] } $e0 = $e1 } $lambda } # 逆ベキ乗法で最小固有値を求める $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) $x = (1.0 ,0.0 ,0.0 ,0.0) # LU分解 forward_elimination $a # 逆ベキ乗法 $lambda = (inverse $a $x) Write-Host Write-Host "eigenvalue" Write-Host ([String]::Format("{0,14:F10}", $lambda)) Write-Host "eigenvector" disp_vector $x
set-variable -name N -value 3 -option constant # 対角要素を表示 function disp_eigenvalue($a) { foreach ($i in 0..$N) { Write-Host ([String]::Format("{0,14:F10}`t", $a[$i][$i])) -nonewline } Write-Host } # LU分解 function decomp($a, $l, $u) { foreach ($i in 0..$N) { foreach ($j in 1..$N) { $l[$i][$j] = 0.0 $u[$i][$j] = 0.0 } } $l[0][0] = 1.0 foreach ($j in 0..$N) { $u[0][$j] = $a[0][$j] } foreach ($i in 1..$N) { $u[$i][0] = 0.0 $l[0][$i] = 0.0 $l[$i][0] = $a[$i][0] / $u[0][0] } foreach ($i in 1..$N) { $l[$i][$i] = 1.0 $t = $a[$i][$i] foreach ($k in 0..$i) { $t -= $l[$i][$k] * $u[$k][$i] } $u[$i][$i] = $t if ($i -lt $N) { foreach ($j in ($i + 1)..$N) { $u[$j][$i] = 0 $l[$i][$j] = 0 $t = $a[$j][$i] foreach ($k in 0..$i) { $t -= $l[$j][$k] * $u[$k][$i] } $l[$j][$i] = $t / $u[$i][$i] $t = $a[$i][$j] foreach ($k in 0..$i) { $t -= $l[$i][$k] * $u[$k][$j] } $u[$i][$j] = $t } } } } # 行列の積 function multiply($a, $b, $c) { foreach ($i in 0..$N) { foreach ($j in 0..$N) { $s = 0.0 foreach ($k in 0..$N) { $s += $a[$i][$k] * $b[$k][$j] } $c[$i][$j] = $s } } } # LR分解で固有値を求める $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) $l = (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0) $u = (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0) foreach ($k in 1..200) { # LU分解 decomp $a $l $u # 行列の積 multiply $u $l $a # 対角要素を表示 Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_eigenvalue $a # 収束判定 $e = 0.0 foreach ($i in 1..$N) { foreach ($j in 0..($i - 1)) { $e += [Math]::Abs($a[$i][$j]) } } if ($e -lt 0.00000000001) { break } } Write-Host Write-Host "固有値" disp_eigenvalue $a
set-variable -name N -value 3 -option constant # 対角要素を表示 function disp_eigenvalue($a) { foreach ($i in 0..$N) { Write-Host ([String]::Format("{0,14:F10}`t", $a[$i][$i])) -nonewline } Write-Host } # QR分解 function decomp($a, $q, $r) { $x = (0.0, 0.0, 0.0, 0.0) foreach ($k in 0..$N) { foreach ($i in 0..$N) { $x[$i] = $a[$i][$k] } if ($k -gt 0) { foreach ($j in 0..($k - 1)) { $t = 0 foreach ($i in 0..$N) { $t += $a[$i][$k] * $q[$i][$j] } $r[$j][$k] = $t $r[$k][$j] = 0 foreach ($i in 0..$N) { $x[$i] -= $t * $q[$i][$j] } } } $t = 0 foreach ($i in 0..$N) { $t += $x[$i] * $x[$i] } $r[$k][$k] = [Math]::Sqrt($t) foreach ($i in 0..$N) { $q[$i][$k] = $x[$i] / $r[$k][$k] } } } # 行列の積 function multiply($a, $b, $c) { foreach ($i in 0..$N) { foreach ($j in 0..$N) { $s = 0.0 foreach ($k in 0..$N) { $s += $a[$i][$k] * $b[$k][$j] } $c[$i][$j] = $s } } } # QR分解で固有値を求める $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) $q = (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0) $r = (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0) foreach ($k in 1..200) { # QR分解 decomp $a $q $r # 行列の積 multiply $r $q $a # 対角要素を表示 Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_eigenvalue $a # 収束判定 $e = 0.0 foreach ($i in 1..$N) { foreach ($j in 0..($i - 1)) { $e += [Math]::Abs($a[$i][$j]) } } if ($e -lt 0.00000000001) { break } } Write-Host Write-Host "固有値" disp_eigenvalue $a
set-variable -name N -value 3 -option constant # 対角要素を表示 function disp_eigenvalue($a) { foreach ($i in 0..$N) { Write-Host ([String]::Format("{0,14:F10}`t", $a[$i][$i])) -nonewline } Write-Host } # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 正規化 (ベクトルの長さを1にする) function normarize($x) { $s = 0.0 foreach ($i in 0..$N) { $s += $x[$i] * $x[$i] } $s = [Math]::Sqrt($s) foreach ($i in 0..$N) { $x[$i] /= $s } } # 固有ベクトルを表示 function disp_eigenvector($matrix) { foreach ($row in $matrix) { normarize $row disp_vector $row } } # ヤコビ法 function jacobi($a, $v) { foreach ($k in 1..100) { # 最大値を探す $max_val = 0.0 foreach ($i in 0..($N - 1)) { foreach ($j in ($i + 1)..$N) { if ($max_val -lt [Math]::Abs($a[$i][$j])) { $max_val = [Math]::Abs($a[$i][$j]) $p = $i $q = $j } } } # θ を求める if ([Math]::Abs($a[$p][$p] - $a[$q][$q]) -lt 0.00000000001) { # a_{pp} = a_{qq} のとき、回転角tをπ/4にする $t = [Math]::PI / 4.0 if ($a[$p][$p] -lt 0.0) { $t = -$t } } else { # a_{pp} ≠ a_{qq} のとき $t = [Math]::Atan(2.0 * $a[$p][$q] / ($a[$p][$p] - $a[$q][$q])) / 2.0 } # θ を使って 行列 U を作成し、A = U^t × A × U $c = [Math]::Cos($t) $s = [Math]::Sin($t) # U^t × A foreach ($i in 0..$N) { $t1 = $a[$p][$i] * $c + $a[$q][$i] * $s $t2 = -$a[$p][$i] * $s + $a[$q][$i] * $c $a[$p][$i] = $t1 $a[$q][$i] = $t2 # 固有ベクトル $t1 = $v[$p][$i] * $c + $v[$q][$i] * $s $t2 = -$v[$p][$i] * $s + $v[$q][$i] * $c $v[$p][$i] = $t1 $v[$q][$i] = $t2 } # A × U foreach ($i in 0..$N) { $t1 = $a[$i][$p] * $c + $a[$i][$q] * $s $t2 = -$a[$i][$p] * $s + $a[$i][$q] * $c $a[$i][$p] = $t1 $a[$i][$q] = $t2 } # 行列の対角要素を表示 (固有値) Write-Host ([String]::Format("{0,3:D}`t", $k)) -nonewline disp_eigenvalue $a # 収束判定 if ($max_val -lt 0.00000000001) { break } } } # ヤコビ法で固有値を求める $a = (5.0, 4.0, 1.0, 1.0), (4.0, 5.0, 1.0, 1.0), (1.0, 1.0, 4.0, 2.0), (1.0, 1.0, 2.0, 4.0) $v = (1.0 ,0.0 ,0.0 ,0.0), (0.0 ,1.0 ,0.0 ,0.0), (0.0 ,0.0 ,1.0 ,0.0), (0.0 ,0.0 ,0.0 ,1.0) # ヤコビ法 jacobi $a $v Write-Host Write-Host "固有値" disp_eigenvalue $a Write-Host Write-Host "固有ベクトル" disp_eigenvector $v
set-variable -name N -value 3 -option constant # 1次元配列を表示 function disp_vector($row) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } # 2次元配列を表示 function disp_matrix($matrix) { foreach ($row in $matrix) { foreach ($col in $row) { Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline } Write-Host } } # ハウスホルダー変換 function tridiagonalize($a, $d, $e) { $v = (0.0 ,0.0 ,0.0 ,0.0) foreach ($k in 0..($N - 2)) { $d[$k] = $a[$k][$k] $t = 0.0 $w = (0.0, 0.0, 0.0, 0.0) foreach ($i in ($k + 1)..$N) { $w[$i] = $a[$i][$k] $t += ($w[$i] * $w[$i]) } $s = [Math]::Sqrt($t) if ($w[$k + 1] -lt 0) { $s = -$s } if ([Math]::Abs($s) -lt 0.00000000001) { $e[$k + 1] = 0.0 } else { $e[$k + 1] = -$s $w[$k + 1] += $s $s = 1 / [Math]::Sqrt($w[$k + 1] * $s) foreach ($i in ($k + 1)..$N) { $w[$i] *= $s } foreach ($i in ($k + 1)..$N) { $s = 0.0 foreach ($j in ($k + 1)..$N) { if ($j -le $i) { $s += $a[$i][$j] * $w[$j] } else { $s += $a[$j][$i] * $w[$j] } } $v[$i] = $s } $s = 0 foreach ($i in ($k + 1)..$N) { $s += $w[$i] * $v[$i] } $s /= 2 foreach ($i in ($k + 1)..$N) { $v[$i] -= $s * $w[$i] } foreach ($i in ($k + 1)..$N) { foreach ($j in ($k + 1)..$i) { $a[$i][$j] -= $w[$i] * $v[$j] + $w[$j] * $v[$i] } } # 次の行は固有ベクトルを求めないなら不要 foreach ($i in ($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] # 次の行は固有ベクトルを求めないなら不要 foreach ($k in ($N..0)) { $w = (0.0, 0.0, 0.0, 0.0) if ($k -lt ($N - 1)) { foreach ($i in ($k + 1)..$N) { $w[$i] = $a[$i][$k] } foreach ($i in ($k + 1)..$N) { $s = 0.0 foreach ($j in ($k + 1)..$N) { $s += $a[$i][$j] * $w[$j] } $v[$i] = $s } foreach ($i in ($k + 1)..$N) { foreach ($j in ($k + 1)..$N) { $a[$i][$j] -= $v[$i] * $w[$j] } } } foreach ($i in 0..$N) { $a[$i][$k] = 0.0 } $a[$k][$k] = 1.0 } } # QR分解 function decomp($a, $d, $e) { $e[0] = 1.0 $h = $N while ([Math]::Abs($e[$h]) -lt 0.00000000001) { $h-- } while ($h -gt 0) { $e[0] = 0.0 $l = $h - 1 while ([Math]::Abs($e[$l]) -ge 0.00000000001) { $l-- } foreach ($j in 1..100) { $w = ($d[$h - 1] - $d[$h]) / 2.0 $s = [Math]::Sqrt($w * $w + $e[$h] * $e[$h]) if ($w -lt 0.0) { $s = - $s } $x = $d[$l] - $d[$h] + $e[$h] * $e[$h] / ($w + $s) $y = $e[$l + 1] $z = 0.0 foreach ($k in $l..($h - 1)) { if ([Math]::Abs($x) -ge [Math]::Abs($y)) { $t = -$y / $x $u = 1 / [Math]::Sqrt($t * $t + 1.0) $s = $t * $u } else { $t = -$x / $y $s = 1 / [Math]::Sqrt($t * $t + 1.0) if ($t -lt 0) { $s = -$s } $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 # 次の行は固有ベクトルを求めないなら不要 foreach ($i in 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 -lt $N - 1) { $x = $e[$k + 1] $y = -$s * $e[$k + 2] $z = $y $e[$k + 2] = $u * $e[$k + 2] } } Write-Host ([String]::Format("{0,3:D}`t", $j)) -nonewline disp_vector $d # 収束判定 if ([Math]::Abs($e[$h]) -lt 0.00000000001) { break } } $e[0] = 1.0 while ([Math]::Abs($e[$h]) -lt 0.00000000001) { $h-- } } # 次の行は固有ベクトルを求めないなら不要 foreach ($k in 0..($N - 1)) { $l = $k foreach ($i in ($k + 1)..$N) { if ($d[$i] -gt $d[$l]) { $l = $i } } $t = $d[$k] $d[$k] = $d[$l] $d[$l] = $t foreach ($i in 0..$N) { $t = $a[$k][$i] $a[$k][$i] = $a[$l][$i] $a[$l][$i] = $t } } } # ハウスホルダー変換とQR分解で固有値を求める $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) $d = (0.0 ,0.0 ,0.0 ,0.0) $e = (0.0 ,0.0 ,0.0 ,0.0) # ハウスホルダー変換 tridiagonalize $a $d $e # QR分解 decomp $a $d $e Write-Host Write-Host "固有値" disp_vector $d Write-Host Write-Host "固有ベクトル" disp_matrix $a