さまざまな言語で数値計算
Only Do What Only You Can Do
ルンゲ・クッタ・ギル法
初期値 $ x_0 $ から次の式によって, 順次 $ x_1, x_2, \dots $ を求める.
例題として, 初速 $ 250 \mathrm{km/h}, 45^\circ $ の角度で打ったボールの軌跡を, ルンゲ・クッタ・ギル法で計算する
(空気抵抗係数を 0.01 で計算)
本来は一時変数をうまく使って, 必要なメモリを少なく抑えるためのアルゴリズムらしいが,
今回はそこまでやらない.
重力による鉛直方向の減速分は, 重力加速度を $g$, 時間を $t$ とすると,
空気抵抗による水平方向の減速分は,速度を $v$, 速度の水平方向成分を $vx$, 空気抵抗係数を $k$ とすると,
同様に, 鉛直方向の減速分は, 速度の鉛直方向成分を $vy$ とすると,
VBScript
Option Explicit Private Const PI = 3.14159265359 '重力加速度 Private Const g = -9.8 '空気抵抗係数 Private Const k = -0.01 '時間間隔(秒) Private Const h = 0.01 '角度 Private Const degree = 45 Private radian: radian = degree * PI / 180.0 '初速 250 km/h -> 秒速に変換 Private v: v = 250 * 1000 \ 3600 '水平方向の速度 Private vx(): ReDim vx(5) vx(0) = v * Cos(radian) '鉛直方向の速度 Private vy(): ReDim vy(5) vy(0) = v * Sin(radian) '経過秒数 Private t: t = 0.0 '位置 Private x(): ReDim x(5) x(0) = 0.0 Private y(): ReDim y(5) y(0) = 0.0 '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx, ByVal vy) fx = k * Sqr(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx, ByVal vy) fy = g + (k * Sqr(vx * vx + vy * vy) * vy) End Function 'Runge-Kutta-Gill法 Dim i: i = 1 Do While (y(0) >= 0.0) '経過秒数 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)) Dim wx: wx = vx(0) + vx(1) / 2.0 Dim wy: wy = vy(0) + vy(1) / 2.0 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) * ((Sqr(2.0) - 1) / 2) + vx(2) * (1 - 1 / Sqr(2.0)) wy = vy(0) + vy(1) * ((Sqr(2.0) - 1) / 2) + vy(2) * (1 - 1 / Sqr(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) / Sqr(2.0) + vx(3) * (1 + 1 / Sqr(2.0)) wy = vy(0) - vy(2) / Sqr(2.0) + vy(3) * (1 + 1 / Sqr(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(0) + ( x(1) + x(2) * (2 - Sqr(2.0)) + x(3) * (2 + Sqr(2.0)) + x(4)) / 6 y(0) = y(0) + ( y(1) + y(2) * (2 - Sqr(2.0)) + y(3) * (2 + Sqr(2.0)) + y(4)) / 6 vx(0) = vx(0) + (vx(1) + vx(2) * (2 - Sqr(2.0)) + vx(3) * (2 + Sqr(2.0)) + vx(4)) / 6 vy(0) = vy(0) + (vy(1) + vy(2) * (2 - Sqr(2.0)) + vy(3) * (2 + Sqr(2.0)) + vy(4)) / 6 WScript.StdOut.Write Right(Space(4) & FormatNumber(t, 2, -1, 0, 0), 4) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(vx(0), 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(vy(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(x(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.WriteLine Right(Space(9) & FormatNumber(y(0), 5, -1, 0, 0), 9) i = i + 1 Loop
Z:\>cscript //nologo Z:\0805.vbs 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
JScript
// 重力加速度 var g = -9.8 // 空気抵抗係数 var k = -0.01 // 時間間隔(秒) var h = 0.01 // 角度 var degree = 45 var radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 var v = parseInt(250 * 1000 / 3600) // 水平方向の速度 var vx = [] vx[0] = v * Math.cos(radian) // 鉛直方向の速度 var vy = [] vy[0] = v * Math.sin(radian) // 経過秒数 var t = 0.0 // 位置 var x = [] x[0] = 0.0 var y = [] y[0] = 0.0 // Runge-Kutta-Gill法 for (var 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]) var wx = vx[0] + vx[1] / 2 var 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 WScript.StdOut.Write((" " + t.toFixed(2) ).slice(-4) + "\t") WScript.StdOut.Write((" " + vx[0].toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + vy[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + x[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + y[0].toFixed(5) ).slice(-8) + "\n") } // 空気抵抗による水平方向の減速分 function fx(vx, vy) { return k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 function fy(vx, vy) { return g + (k * Math.sqrt(vx * vx + vy * vy) * vy) }
Z:\>cscript //nologo Z:\0805.js 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
PowerShell
# 重力加速度 $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])) }
Z:\>powershell -file Z:\0805.ps1 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Perl
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); }
Z:\>perl Z:\0805.pl 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
PHP
<?php # 重力加速度 $g = -9.8; # 空気抵抗係数 $k = -0.01; # 時間間隔(秒) $h = 0.01; # 角度 $degree = 45; $radian = $degree * M_PI / 180.0; # 初速 250 km/h -> 秒速に変換 $v = (int)(250 * 1000 / 3600); # 水平方向の速度 $vx = array(); $vx[0] = $v * cos($radian); # 鉛直方向の速度 $vy = array(); $vy[0] = $v * sin($radian); # 経過秒数 $t = 0.0; # 位置 $x = array(); $x[0] = 0.0; $y = array(); $y[0] = 0.0; # Runge-Kutta-Gill法 for ($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]); $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] * ((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]); } # 空気抵抗による水平方向の減速分 function fx($vx, $vy) { global $k; return $k * sqrt($vx * $vx + $vy * $vy) * $vx; } # 重力と空気抵抗による鉛直方向の減速分 function fy($vx, $vy) { global $g, $k; return $g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy); } ?>
Z:\>php Z:\0805.php 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Python
# coding: Shift_JIS import math # 重力加速度 g = -9.8 # 空気抵抗係数 k = -0.01 # 時間間隔(秒) h = 0.01 # 角度 degree = 45 radian = degree * math.pi / 180.0 # 初速 250 km/h -> 秒速に変換 v = int(250 * 1000 / 3600) # 水平方向の速度 vx = [0 for i in range(5)] vx[0] = v * math.cos(radian) # 鉛直方向の速度 vy = [0 for i in range(5)] vy[0] = v * math.sin(radian) # 経過秒数 t = 0.0 # 位置 x = [0 for i in range(5)] x[0] = 0.0 y = [0 for i in range(5)] y[0] = 0.0 # 空気抵抗による水平方向の減速分 def fx(vx, vy): return k * math.sqrt(vx * vx + vy * vy) * vx # 重力と空気抵抗による鉛直方向の減速分 def fy(vx, vy): return g + (k * math.sqrt(vx * vx + vy * vy) * vy) # Runge-Kutta-Gill法 i = 1 while y[0] >= 0.0: # 経過秒数 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 print "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" % (t, vx[0], vy[0], x[0], y[0]) i += 1
Z:\>python Z:\0805.py 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Ruby
# 重力加速度 $g = -9.8 # 空気抵抗係数 $k = -0.01 # 時間間隔(秒) h = 0.01 # 角度 degree = 45 radian = degree * Math::PI / 180.0 # 初速 250 km/h -> 秒速に変換 v = 250 * 1000 / 3600 # 水平方向の速度 vx = Array.new(5) vx[0] = v * Math.cos(radian) # 鉛直方向の速度 vy = Array.new(5) vy[0] = v * Math.sin(radian) # 経過秒数 t = 0.0 # 位置 x = Array.new(5) x[0] = 0.0 y = Array.new(5) y[0] = 0.0 # 空気抵抗による水平方向の減速分 def fx(vx, vy) return $k * Math.sqrt(vx * vx + vy * vy) * vx end # 重力と空気抵抗による鉛直方向の減速分 def fy(vx, vy) return $g + ($k * Math.sqrt(vx * vx + vy * vy) * vy) end # Runge-Kutta-Gill法 i = 1 while y[0] >= 0.0 do # 経過秒数 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 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]) i += 1 end
Z:\>ruby Z:\0805.rb 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Groovy
Pascal
program Pas0805(arg); {$MODE delphi} uses SysUtils, Math; const // 重力加速度 g = -9.8; // 空気抵抗係数 k = -0.01; // 時間間隔(秒) h = 0.01; // 角度 degree = 45; // 空気抵抗による水平方向の減速分 function fx(vx:Double; vy:Double):Double; begin result := k * Sqrt(vx * vx + vy * vy) * vx; end; // 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Double; vy:Double):Double; begin result := g + (k * Sqrt(vx * vx + vy * vy) * vy); end; var // 角度 radian:Double; // 初速 v:Double; // 水平方向の速度 vx:array [0..4] of Double; wx:Double; // 鉛直方向の速度 vy:array [0..4] of Double; wy:Double; // 経過秒数 t:Double = 0.0; // 位置 x:array [0..4] of Double; y:array [0..4] of Double; i:Integer; begin // 角度 radian := degree * PI / 180.0; // 初速 250 km/h -> 秒速に変換 v := 250 * 1000 div 3600; // 水平方向の速度 vx[0] := v * Cos(radian); // 鉛直方向の速度 vy[0] := v * Sin(radian); // 位置 x[0] := 0.0; y[0] := 0.0; // Runge-Kutta-Gill法 i := 1; while y[0] >= 0.0 do begin // 経過秒数 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] * ((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[0] + ( x[1] + x[2] * (2 - Sqrt(2.0)) + x[3] * (2 + Sqrt(2.0)) + x[4]) / 6; y[0] := y[0] + ( y[1] + y[2] * (2 - Sqrt(2.0)) + y[3] * (2 + Sqrt(2.0)) + y[4]) / 6; vx[0] := vx[0] + (vx[1] + vx[2] * (2 - Sqrt(2.0)) + vx[3] * (2 + Sqrt(2.0)) + vx[4]) / 6; vy[0] := vy[0] + (vy[1] + vy[2] * (2 - Sqrt(2.0)) + vy[3] * (2 + Sqrt(2.0)) + vy[4]) / 6; writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%9.5f', [t, vx[0], vy[0], x[0], y[0]])); inc(i); end; end.
Z:\>fpc -v0 -l- Pas0805.pp Z:\>Pas0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Ada
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; procedure Ada0805 is -- 重力加速度 g : Constant Long_Float := -9.8; -- 空気抵抗係数 k : Constant Long_Float := -0.01; -- 時間間隔(秒) h : Constant Long_Float := 0.01; -- 角度 degree : Constant Long_Float := 45.0; -- 空気抵抗による水平方向の減速分 function fx(vx:Long_Float; vy:Long_Float) return Long_Float is begin return k * Sqrt(vx * vx + vy * vy) * vx; end fx; -- 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Long_Float; vy:Long_Float) return Long_Float is begin return g + (k * Sqrt(vx * vx + vy * vy) * vy); end fy; -- 角度 radian:Long_Float; -- 初速 v:Long_Float; -- 水平方向の速度 vx:array (0..4) of Long_Float; wx:Long_Float; -- 鉛直方向の速度 vy:array (0..4) of Long_Float; wy:Long_Float; -- 経過秒数 t:Long_Float := 0.0; -- 位置 x:array (0..4) of Long_Float; y:array (0..4) of Long_Float; i:Integer; begin -- 角度 radian := degree * Pi / 180.0; -- 初速 250 km/h -> 秒速に変換 v := Long_Float(250 * 1000 / 3600); -- 水平方向の速度 vx(0) := v * Cos(radian); -- 鉛直方向の速度 vy(0) := v * Sin(radian); -- 位置 x(0) := 0.0; y(0) := 0.0; -- Runge-Kutta-Gill法 i := 1; while y(0) >= 0.0 loop -- 経過秒数 t := Long_Float(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.0; wy := vy(0) + vy(1) / 2.0; 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.0) / 2.0) + vx(2) * (1.0 - 1.0 / Sqrt(2.0)); wy := vy(0) + vy(1) * ((Sqrt(2.0) - 1.0) / 2.0) + vy(2) * (1.0 - 1.0 / 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.0 + 1.0 / Sqrt(2.0)); wy := vy(0) - vy(2) / Sqrt(2.0) + vy(3) * (1.0 + 1.0 / 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(0) + ( x(1) + x(2) * (2.0 - Sqrt(2.0)) + x(3) * (2.0 + Sqrt(2.0)) + x(4)) / 6.0; y(0) := y(0) + ( y(1) + y(2) * (2.0 - Sqrt(2.0)) + y(3) * (2.0 + Sqrt(2.0)) + y(4)) / 6.0; vx(0) := vx(0) + (vx(1) + vx(2) * (2.0 - Sqrt(2.0)) + vx(3) * (2.0 + Sqrt(2.0)) + vx(4)) / 6.0; vy(0) := vy(0) + (vy(1) + vy(2) * (2.0 - Sqrt(2.0)) + vy(3) * (2.0 + Sqrt(2.0)) + vy(4)) / 6.0; Put(t, Fore=>1, Aft=>2, Exp=>0); Put(Ascii.HT); Put(vx(0), Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(vy(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(x(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(y(0), Fore=>4, Aft=>5, Exp=>0); New_Line; i := i + 1; end loop; end Ada0805;
xxxxxx@yyyyyy /Z $ gnatmake Ada0805.adb xxxxxx@yyyyyy /Z $ Ada0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
VB.NET
Option Explicit Module VB0805 '重力加速度 Private Const g As Double = -9.8 '空気抵抗係数 Private Const k As Double = -0.01 '時間間隔(秒) Private Const h As Double = 0.01 Public Sub Main() '角度 Const degree As Double = 45 Dim radian As Double = degree * Math.PI / 180.0 '初速 250 km/h -> 秒速に変換 Dim v As Double = 250 * 1000 \ 3600 '水平方向の速度 Dim vx(5) As Double vx(0) = v * Math.Cos(radian) '鉛直方向の速度 Dim vy(5) As Double vy(0) = v * Math.Sin(radian) '経過秒数 Dim t As Double = 0.0 '位置 Dim x(5) As Double x(0) = 0.0 Dim y(5) As Double y(0) = 0.0 'Runge-Kutta-Gill法 Dim i As Integer = 1 Do While (y(0) >= 0.0) '経過秒数 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)) Dim wx As Double = vx(0) + vx(1) / 2.0 Dim wy As Double = vy(0) + vy(1) / 2.0 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 Console.WriteLine(String.Format("{0,4:F2}{5}{1,8:F5}{5}{2,9:F5}{5}{3,9:F5}{5}{4,9:F5}", t, vx(0), vy(0), x(0), y(0), vbTab)) i += 1 Loop End Sub '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx As Double, ByVal vy As Double) As Double Return k * Math.Sqrt(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx As Double, ByVal vy As Double) As Double Return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) End Function End Module
Z:\>vbc -nologo VB0805.vb Z:\>VB0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
C#
using System; public class CS0805 { // 重力加速度 private const double g = -9.8; // 空気抵抗係数 private const double k = -0.01; // 時間間隔(秒) private const double h = 0.01; public static void Main() { // 角度 double degree = 45; double radian = degree * Math.PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double[] vx = new double[5]; vx[0] = v * Math.Cos(radian); // 鉛直方向の速度 double[] vy = new double[5]; vy[0] = v * Math.Sin(radian); // 経過秒数 double t = 0.0; // 位置 double[] x = new double[5]; x[0] = 0.0; double[] y = new double[5]; y[0] = 0.0; // Runge-Kutta-Gill法 for (int 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]); double wx = vx[0] + vx[1] / 2; double 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; Console.WriteLine(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])); } } // 空気抵抗による水平方向の減速分 private static double fx(double vx, double vy) { return k * Math.Sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 private static double fy(double vx, double vy) { return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy); } }
Z:\>csc -nologo CS0805.cs Z:\>CS0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Java
import static java.lang.System.out; public class Java0805 { // 重力加速度 private static final double g = -9.8; // 空気抵抗係数 private static final double k = -0.01; // 時間間隔(秒) private static final double h = 0.01; public static void main(String []args) { // 角度 double degree = 45; double radian = degree * Math.PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double[] vx = new double[5]; vx[0] = v * Math.cos(radian); // 鉛直方向の速度 double[] vy = new double[5]; vy[0] = v * Math.sin(radian); // 経過秒数 double t = 0.0; // 位置 double[] x = new double[5]; x[0] = 0.0; double[] y = new double[5]; y[0] = 0.0; // Runge-Kutta-Gill法 for (int 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]); double wx = vx[0] + vx[1] / 2; double 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; out.println(String.format("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f", t, vx[0], vy[0], x[0], y[0])); } } // 空気抵抗による水平方向の減速分 private static double fx(double vx, double vy) { return k * Math.sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 private static double fy(double vx, double vy) { return g + (k * Math.sqrt(vx * vx + vy * vy) * vy); } }
Z:\>javac Java0805.java Z:\>java Java0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
C++
#include <iostream> #include <iomanip> #include <math.h> using namespace std; // 重力加速度 const double g = -9.8; // 空気抵抗係数 const double k = -0.01; // 時間間隔(秒) const double h = 0.01; // 空気抵抗による水平方向の減速分 double fx(double vx, double vy); // 重力と空気抵抗による鉛直方向の減速分 double fy(double vx, double vy); int main() { // 角度 double degree = 45; double radian = degree * M_PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double vx[5]; vx[0] = v * cos(radian); // 鉛直方向の速度 double vy[5]; vy[0] = v * sin(radian); // 経過秒数 double t = 0.0; // 位置 double x[5]; x[0] = 0.0; double y[5]; y[0] = 0.0; // Runge-Kutta-Gill法 for (int i = 1; y[0] >= 0.0; i++) { // 経過秒数 t = i * h; cout << setw(4) << fixed << setprecision(2) << t << "\t"; // 速度 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]); double wx = vx[0] + vx[1] / 2; double 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; cout << setw(8) << fixed << setprecision(5) << vx[0] << "\t"; cout << setw(9) << fixed << setprecision(5) << vy[0] << "\t"; cout << setw(9) << fixed << setprecision(5) << x[0] << "\t"; cout << setw(8) << fixed << setprecision(5) << y[0] << endl; } return 0; } // 空気抵抗による水平方向の減速分 double fx(double vx, double vy) { return k * sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 double fy(double vx, double vy) { return g + (k * sqrt(vx * vx + vy * vy) * vy); }
Z:\>bcc32 -q CP0805.cpp cp0805.cpp: Z:\>CP0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Objective-C
#import <Foundation/Foundation.h> #import <math.h> // 重力加速度 const double g = -9.8; // 空気抵抗係数 const double k = -0.01; // 時間間隔(秒) const double h = 0.01; // 空気抵抗による水平方向の減速分 double fx(double vx, double vy); // 重力と空気抵抗による鉛直方向の減速分 double fy(double vx, double vy); int main() { // 角度 double degree = 45; double radian = degree * M_PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double vx[5]; vx[0] = v * cos(radian); // 鉛直方向の速度 double vy[5]; vy[0] = v * sin(radian); // 経過秒数 double t = 0.0; // 位置 double x[5]; x[0] = 0.0; double y[5]; y[0] = 0.0; // Runge-Kutta-Gill法 int i; for (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]); double wx = vx[0] + vx[1] / 2; double 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]); } return 0; } // 空気抵抗による水平方向の減速分 double fx(double vx, double vy) { return k * sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 double fy(double vx, double vy) { return g + (k * sqrt(vx * vx + vy * vy) * vy); }
xxxxxx@yyyyyy /Z $ gcc -o OC0805 OC0805.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
D
import std.stdio; import std.math; // 重力加速度 const double g = -9.8; // 空気抵抗係数 const double k = -0.01; // 時間間隔(秒) const double h = 0.01; void main(string[] args) { // 角度 double degree = 45; double radian = degree * PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double vx[5]; vx[0] = v * cos(radian); // 鉛直方向の速度 double vy[5]; vy[0] = v * sin(radian); // 経過秒数 double t = 0.0; // 位置 double x[5]; x[0] = 0.0; double y[5]; y[0] = 0.0; // Runge-Kutta-Gill法 for (int 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]); double wx = vx[0] + vx[1] / 2; double 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; writefln("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f", t, vx[0], vy[0], x[0], y[0]); } } // 空気抵抗による水平方向の減速分 double fx(double vx, double vy) { return k * sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 double fy(double vx, double vy) { return g + (k * sqrt(vx * vx + vy * vy) * vy); }
Z:\>dmd D0805.d Z:\>D0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02079 6.31 8.99347 -24.15100 126.75335 -0.22053
Go
package main import "fmt" import "math" // 重力加速度 const g float64 = -9.8 // 空気抵抗係数 const k float64 = -0.01 // 時間間隔(秒) const h float64 = 0.01 func main() { // 角度 var degree float64 = 45 var radian float64 = degree * math.Pi / 180.0 // 初速 250 km/h -> 秒速に変換 var v float64 = 250 * 1000 / 3600 // 水平方向の速度 var vx[5] float64 vx[0] = v * math.Cos(radian) // 鉛直方向の速度 var vy[5] float64 vy[0] = v * math.Sin(radian) // 経過秒数 var t float64 = 0.0 // 位置 var x[5] float64 x[0] = 0.0 var y[5] float64 y[0] = 0.0 // Runge-Kutta-Gill法 for i := 1; y[0] >= 0.0; i++ { // 経過秒数 t = float64(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]) var wx float64 = vx[0] + vx[1] / 2 var wy float64 = 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 fmt.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]) } } // 空気抵抗による水平方向の減速分 func fx(vx float64, vy float64) float64 { return k * math.Sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 func fy(vx float64, vy float64) float64 { return g + (k * math.Sqrt(vx * vx + vy * vy) * vy) }
Z:\>8g GO0805.go Z:\>8l -o GO0805.exe GO0805.8 Z:\>GO0805 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Scala
object Scala0804 { // 重力加速度 val g = -9.8 // 空気抵抗係数 val k = -0.01 // 時間間隔(秒) val h = 0.01 def main(args: Array[String]) { // 角度 val degree = 45 val radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 val v = 250 * 1000 / 3600 // 水平方向の速度 val vx = v * Math.cos(radian) // 鉛直方向の速度 val vy = v * Math.sin(radian) // 位置 val x = 0.0 val y = 0.0 // Runge-Kutta-Gill法 rungekuttagill(1, vx, vy, x, y) } def rungekuttagill(i:Int, vx:Double, vy:Double, x:Double, y:Double):Unit = { // 経過秒数 val t = i * h // 位置・速度 val wx1 = h * vx val wy1 = h * vy val wvx1 = h * fx(vx, vy) val wvy1 = h * fy(vx, vy) val wvx5 = vx + wvx1 / 2 val wvy5 = vy + wvy1 / 2 val wx2 = h * wvx5 val wy2 = h * wvy5 val wvx2 = h * fx(wvx5, wvy5) val wvy2 = h * fy(wvx5, wvy5) val wvx6 = vx + wvx1 * ((Math.sqrt(2.0) - 1) / 2) + wvx2 * (1 - 1 / Math.sqrt(2.0)) val wvy6 = vy + wvy1 * ((Math.sqrt(2.0) - 1) / 2) + wvy2 * (1 - 1 / Math.sqrt(2.0)) val wx3 = h * wvx6 val wy3 = h * wvy6 val wvx3 = h * fx(wvx6, wvy6) val wvy3 = h * fy(wvx6, wvy6) val wvx7 = vx - wvx2 / Math.sqrt(2.0) + wvx3 * (1 + 1 / Math.sqrt(2.0)) val wvy7 = vy - wvy2 / Math.sqrt(2.0) + wvy3 * (1 + 1 / Math.sqrt(2.0)) val wx4 = h * wvx7 val wy4 = h * wvy7 val wvx4 = h * fx(wvx7, wvy7) val wvy4 = h * fy(wvx7, wvy7) val wx = x + ( wx1 + wx2 * (2 - Math.sqrt(2.0)) + wx3 * (2 + Math.sqrt(2.0)) + wx4) / 6 val wy = y + ( wy1 + wy2 * (2 - Math.sqrt(2.0)) + wy3 * (2 + Math.sqrt(2.0)) + wy4) / 6 val wvx = vx + (wvx1 + wvx2 * (2 - Math.sqrt(2.0)) + wvx3 * (2 + Math.sqrt(2.0)) + wvx4) / 6 val wvy = vy + (wvy1 + wvy2 * (2 - Math.sqrt(2.0)) + wvy3 * (2 + Math.sqrt(2.0)) + wvy4) / 6 println("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f".format(t, wvx, wvy, wx, wy)) if (wy >= 0.0) rungekuttagill(i+1, wvx, wvy, wx, wy) else () } // 空気抵抗による水平方向の減速分 def fx(vx:Double, vy:Double) = { k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 def fy(vx:Double, vy:Double) = { g + (k * Math.sqrt(vx * vx + vy * vy) * vy) } }
Z:\>scala Scala0805.scala 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
F#
module Fs0805 open System // 重力加速度 let g = -9.8 // 空気抵抗係数 let k = -0.01 // 時間間隔(秒) let h = 0.01 // 角度 let degree = 45.0 let radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 let v = float(250 * 1000 / 3600) // 水平方向の速度 let vx = v * Math.Cos(radian) // 鉛直方向の速度 let vy = v * Math.Sin(radian) // 位置 let x = 0.0 let y = 0.0 // 空気抵抗による水平方向の減速分 let fx(vx:Double) (vy:Double) = k * Math.Sqrt(vx * vx + vy * vy) * vx // 重力と空気抵抗による鉛直方向の減速分 let fy(vx:Double) (vy:Double) = g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) // Runge-Kutta-Gill法 let rec rungekuttagill(i:int) (vx:double) (vy:double) (x:double) (y:double):unit = // 経過秒数 let t = float(i) * h // 位置・速度 let wx1 = h * vx let wy1 = h * vy let wvx1 = h * (fx vx vy) let wvy1 = h * (fy vx vy) let wvx5 = vx + wvx1 / 2.0 let wvy5 = vy + wvy1 / 2.0 let wx2 = h * wvx5 let wy2 = h * wvy5 let wvx2 = h * (fx wvx5 wvy5) let wvy2 = h * (fy wvx5 wvy5) let wvx6 = vx + wvx1 * ((Math.Sqrt(2.0) - 1.0) / 2.0) + wvx2 * (1.0 - 1.0 / Math.Sqrt(2.0)) let wvy6 = vy + wvy1 * ((Math.Sqrt(2.0) - 1.0) / 2.0) + wvy2 * (1.0 - 1.0 / Math.Sqrt(2.0)) let wx3 = h * wvx6 let wy3 = h * wvy6 let wvx3 = h * (fx wvx6 wvy6) let wvy3 = h * (fy wvx6 wvy6) let wvx7 = vx - wvx2 / Math.Sqrt(2.0) + wvx3 * (1.0 + 1.0 / Math.Sqrt(2.0)) let wvy7 = vy - wvy2 / Math.Sqrt(2.0) + wvy3 * (1.0 + 1.0 / Math.Sqrt(2.0)) let wx4 = h * wvx7 let wy4 = h * wvy7 let wvx4 = h * (fx wvx7 wvy7) let wvy4 = h * (fy wvx7 wvy7) let wx = x + ( wx1 + wx2 * (2.0 - Math.Sqrt(2.0)) + wx3 * (2.0 + Math.Sqrt(2.0)) + wx4) / 6.0 let wy = y + ( wy1 + wy2 * (2.0 - Math.Sqrt(2.0)) + wy3 * (2.0 + Math.Sqrt(2.0)) + wy4) / 6.0 let wvx = vx + (wvx1 + wvx2 * (2.0 - Math.Sqrt(2.0)) + wvx3 * (2.0 + Math.Sqrt(2.0)) + wvx4) / 6.0 let wvy = vy + (wvy1 + wvy2 * (2.0 - Math.Sqrt(2.0)) + wvy3 * (2.0 + Math.Sqrt(2.0)) + wvy4) / 6.0 printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy if wy >= 0.0 then (rungekuttagill (i+1) wvx wvy wx wy) else () // Runge-Kutta-Gill法 (rungekuttagill 1 vx vy x y) exit 0
Z:\>fsi --nologo --quiet Fs0805.fs 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Clojure
; 重力加速度 (def g -9.8) ; 空気抵抗係数 (def k -0.01) ; 時間間隔(秒) (def h 0.01) ; 角度 (def degree 45.0) (def radian (* degree (/ (. Math PI) 180.0))) ; 初速 250 km/h -> 秒速に変換 (def v (quot (* 250 1000) 3600)) ; 水平方向の速度 (def vx (* v (. Math cos radian))) ; 鉛直方向の速度 (def vy (* v (. Math sin radian))) ; 位置 (def x 0.0) (def y 0.0) ; 空気抵抗による水平方向の減速分 (defn fx[vx vy] (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx))) ; 重力と空気抵抗による鉛直方向の減速分 (defn fy[vx vy] (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy)))) ; Runge-Kutta-Gill法 (defn rungekuttagill[i vx vy x y] ; 経過秒数 (def t (* i h)) ; 位置・速度 (def wx1 (* h vx)) (def wy1 (* h vy)) (def wvx1 (* h (fx vx vy))) (def wvy1 (* h (fy vx vy))) (def wvx5 (+ vx (/ wvx1 2.0))) (def wvy5 (+ vy (/ wvy1 2.0))) (def wx2 (* h wvx5)) (def wy2 (* h wvy5)) (def wvx2 (* h (fx wvx5 wvy5))) (def wvy2 (* h (fy wvx5 wvy5))) (def wvx6 (+ vx (* wvx1 (/ (- (. Math sqrt 2.0) 1.0) 2.0)) (* wvx2 (- 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wvy6 (+ vy (* wvy1 (/ (- (. Math sqrt 2.0) 1.0) 2.0)) (* wvy2 (- 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wx3 (* h wvx6)) (def wy3 (* h wvy6)) (def wvx3 (* h (fx wvx6 wvy6))) (def wvy3 (* h (fy wvx6 wvy6))) (def wvx7 (+ (- vx (/ wvx2 (. Math sqrt 2.0))) (* wvx3 (+ 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wvy7 (+ (- vy (/ wvy2 (. Math sqrt 2.0))) (* wvy3 (+ 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wx4 (* h wvx7)) (def wy4 (* h wvy7)) (def wvx4 (* h (fx wvx7 wvy7))) (def wvy4 (* h (fy wvx7 wvy7))) (def wx (+ x (/ (+ wx1 (* wx2 (- 2.0 (. Math sqrt 2.0))) (* wx3 (+ 2.0 (. Math sqrt 2.0))) wx4) 6.0))) (def wy (+ y (/ (+ wy1 (* wy2 (- 2.0 (. Math sqrt 2.0))) (* wy3 (+ 2.0 (. Math sqrt 2.0))) wy4) 6.0))) (def wvx (+ vx (/ (+ wvx1 (* wvx2 (- 2.0 (. Math sqrt 2.0))) (* wvx3 (+ 2.0 (. Math sqrt 2.0))) wvx4) 6.0))) (def wvy (+ vy (/ (+ wvy1 (* wvy2 (- 2.0 (. Math sqrt 2.0))) (* wvy3 (+ 2.0 (. Math sqrt 2.0))) wvy4) 6.0))) (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy)) (if (>= wy 0.0) (rungekuttagill (+ i 1) wvx wvy wx wy) nil)) (rungekuttagill 1 vx vy x y)
Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0805.clj 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053
Haskell
import Text.Printf -- 重力加速度 g = -9.8 :: Double -- 空気抵抗係数 k = -0.01 :: Double -- 時間間隔(秒) h = 0.01 :: Double -- 空気抵抗による水平方向の減速分 fx::Double->Double->Double fx vx vy = let v = sqrt(vx * vx + vy * vy) in k * v * vx -- 空気抵抗による鉛直方向の減速分 fy::Double->Double->Double fy vx vy = let v = sqrt(vx * vx + vy * vy) in g + (k * v * vy) -- Runge-Kutta-Gill法 rungekuttagill::Integer->Double->Double->Double->Double->IO () rungekuttagill i vx vy x y = let -- 経過秒数 t = (fromIntegral i) * h -- 位置・速度 wx1 = h * vx wy1 = h * vy wvx1 = h * (fx vx vy) wvy1 = h * (fy vx vy) wvx5 = vx + wvx1 / 2.0 wvy5 = vy + wvy1 / 2.0 wx2 = h * wvx5 wy2 = h * wvy5 wvx2 = h * (fx wvx5 wvy5) wvy2 = h * (fy wvx5 wvy5) wvx6 = vx + wvx1 * ((sqrt(2.0) - 1.0) / 2.0) + wvx2 * (1.0 - 1.0 / sqrt(2.0)) wvy6 = vy + wvy1 * ((sqrt(2.0) - 1.0) / 2.0) + wvy2 * (1.0 - 1.0 / sqrt(2.0)) wx3 = h * wvx6 wy3 = h * wvy6 wvx3 = h * (fx wvx6 wvy6) wvy3 = h * (fy wvx6 wvy6) wvx7 = vx - wvx2 / sqrt(2.0) + wvx3 * (1.0 + 1.0 / sqrt(2.0)) wvy7 = vy - wvy2 / sqrt(2.0) + wvy3 * (1.0 + 1.0 / sqrt(2.0)) wx4 = h * wvx7 wy4 = h * wvy7 wvx4 = h * (fx wvx7 wvy7) wvy4 = h * (fy wvx7 wvy7) wx = x + ( wx1 + wx2 * (2.0 - sqrt(2.0)) + wx3 * (2.0 + sqrt(2.0)) + wx4) / 6.0 wy = y + ( wy1 + wy2 * (2.0 - sqrt(2.0)) + wy3 * (2.0 + sqrt(2.0)) + wy4) / 6.0 wvx = vx + (wvx1 + wvx2 * (2.0 - sqrt(2.0)) + wvx3 * (2.0 + sqrt(2.0)) + wvx4) / 6.0 wvy = vy + (wvy1 + wvy2 * (2.0 - sqrt(2.0)) + wvy3 * (2.0 + sqrt(2.0)) + wvy4) / 6.0 in if y >= 0.0 then do printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy rungekuttagill (i+1) wvx wvy wx wy else return () main = do -- 角度 let degree = 45.0 :: Double let radian = degree * pi / 180.0 -- 初速 250 km/h -> 秒速に変換 let v = (fromIntegral (250 * 1000 `div` 3600)) -- 水平方向の速度 let vx = v * cos(radian) -- 鉛直方向の速度 let vy = v * sin(radian) -- 位置 let x = 0.0 let y = 0.0 -- Runge-Kutta-Gill法 rungekuttagill 1 vx vy x y
Z:\>runghc Hs0805.hs 0.01 48.45619 48.35852 0.48623 0.48574 0.02 48.12689 47.93222 0.96914 0.96719 0.03 47.80236 47.51133 1.44878 1.44440 0.04 47.48250 47.09575 1.92520 1.91743 0.05 47.16722 46.68537 2.39845 2.38633 0.06 46.85642 46.28007 2.86856 2.85116 0.07 46.55001 45.87974 3.33559 3.31195 0.08 46.24790 45.48430 3.79958 3.76877 0.09 45.94999 45.09363 4.26056 4.22165 省略 6.20 9.25063 -23.74826 125.74997 2.41410 6.21 9.22707 -23.78565 125.84236 2.17643 6.22 9.20355 -23.82289 125.93451 1.93839 6.23 9.18006 -23.85997 126.02643 1.69997 6.24 9.15661 -23.89690 126.11811 1.46119 6.25 9.13319 -23.93366 126.20956 1.22203 6.26 9.10982 -23.97027 126.30078 0.98251 6.27 9.08647 -24.00673 126.39176 0.74263 6.28 9.06317 -24.04303 126.48251 0.50238 6.29 9.03990 -24.07917 126.57302 0.26177 6.30 9.01667 -24.11516 126.66330 0.02080 6.31 8.99347 -24.15100 126.75335 -0.22053