さまざまな言語で数値計算

オイラー法

初期値 $ 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(1)
vx(0) = v * Cos(radian)
'鉛直方向の速度
Private vy(): ReDim vy(1)
vy(0) = v * Sin(radian)
'経過秒数
Private t: t = 0.0
'位置
Private x: x = 0.0
Private y: y = 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

'Euler法
Dim i: i = 1
Do While (y >= 0.0)
    '経過秒数
    t = i * h

    '位置
    x = x + h * vx(0)
    y = y + h * vy(0)

    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,     5, -1, 0, 0), 9) & vbTab
    WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(y,     5, -1, 0, 0), 8)

    '速度
    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)

    i = i + 1
Loop

Z:\>cscript //nologo Z:\0801.vbs
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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 = 0.0
var y = 0.0

// Euler法
for (var i = 1; y >= 0.0; i++)
{
    // 経過秒数
    t = i * h

    // 位置
    x += h * vx[0]
    y += h * vy[0]

    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.toFixed(5)       ).slice(-9) + "\t")
    WScript.StdOut.Write(("        "  + y.toFixed(5)       ).slice(-8) + "\n")

    // 速度
    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]
}

// 空気抵抗による水平方向の減速分
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:\0801.js
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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[] 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]
}

Z:\>powershell -file Z:\0801.ps1
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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 = 0.0;
my $y = 0.0;

# Euler法
for (my $i = 1; $y >= 0.0; $i++)
{
    # 経過秒数
    $t = $i * $h;

    # 位置
    $x += $h * $vx[0];
    $y += $h * $vy[0];

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", $t, $vx[0], $vy[0], $x, $y);

    # 速度
    $vx[1] = $vx[0] + $h * fx($vx[0], $vy[0]);
    $vy[1] = $vy[0] + $h * fy($vx[0], $vy[0]);
    $vx[0] = $vx[1];
    $vy[0] = $vy[1];
}

# 空気抵抗による水平方向の減速分
sub fx
{
    my ($vx, $vy) = @_;
    $k * sqrt($vx * $vx + $vy * $vy) * $vx;
}
# 重力と空気抵抗による鉛直方向の減速分
sub fy
{
    my ($vx, $vy) = @_;
    $g + ($k * sqrt($vx * $vx + $vy * $vy) * $vy);
}

Z:\>perl Z:\0801.pl
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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 = 0.0;
$y = 0.0;

# Euler法
for ($i = 1; $y >= 0.0; $i++)
{
    # 経過秒数
    $t = $i * $h;

    # 位置
    $x += $h * $vx[0];
    $y += $h * $vy[0];

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", $t, $vx[0], $vy[0], $x, $y);

    # 速度
    $vx[1] = $vx[0] + $h * fx($vx[0], $vy[0]);
    $vy[1] = $vy[0] + $h * fy($vx[0], $vy[0]);
    $vx[0] = $vx[1];
    $vy[0] = $vy[1];
}

# 空気抵抗による水平方向の減速分
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:\0801.php
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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(2)]
vx[0] = v * math.cos(radian)
# 鉛直方向の速度
vy = [0 for i in range(2)]
vy[0] = v * math.sin(radian)
# 経過秒数
t = 0.0
# 位置
x = 0.0
y = 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)

# Euler法
i = 1
while y >= 0.0:
    # 経過秒数
    t = i * h

    # 位置
    x += h * vx[0]
    y += h * vy[0]

    print "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" % (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]

    i += 1

Z:\>python Z:\0801.py
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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(2)
vx[0] = v * Math.cos(radian)
# 鉛直方向の速度
vy = Array.new(2)
vy[0] = v * Math.sin(radian)
# 経過秒数
t = 0.0
# 位置
x = 0.0
y = 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

# Euler法
i = 1
while y >= 0.0 do
    # 経過秒数
    t = i * h

    # 位置
    x += h * vx[0]
    y += h * vy[0]

    printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", t, vx[0], vy[0], x, y)

    # 速度
    vx[1] = vx[0] + h * fx(vx[0], vy[0])
    vy[1] = vy[0] + h * fy(vx[0], vy[0])
    vx[0] = vx[1]
    vy[0] = vy[1]

    i += 1
end

Z:\>ruby Z:\0801.rb
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

Groovy

Pascal

program Pas0801(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..1] of Double;
    // 鉛直方向の速度
    vy:array [0..1] of Double;
    // 経過秒数
    t:Double = 0.0;
    // 位置
    x:Double = 0.0;
    y:Double = 0.0;

    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);

    // Euler法
    i := 1;
    while y >= 0.0 do
    begin
        // 経過秒数
        t := i * h;

        // 位置
        x := x + h * vx[0];
        y := y + h * vy[0];

        writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%8.5f', [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];

        inc(i);
    end;
end.

Z:\>fpc -v0 -l- Pas0801.pp

Z:\>Pas0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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 Ada0801 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..1) of Long_Float;
    -- 鉛直方向の速度
    vy:array (0..1) of Long_Float;
    -- 経過秒数
    t:Long_Float := 0.0;
    -- 位置
    x:Long_Float := 0.0;
    y:Long_Float := 0.0;

    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);

    -- Euler法
    i := 1;
    while y >= 0.0 loop
        -- 経過秒数
        t := Long_Float(i) * h;

        -- 位置
        x := x + h * vx(0);
        y := y + h * vy(0);

        Put(t,       Fore=>1, Aft=>2, Exp=>0);
        Put(Ascii.HT);
        Put(vx(1),   Fore=>3, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(vy(1),   Fore=>4, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(x,       Fore=>4, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(y,       Fore=>3, Aft=>5, Exp=>0);
        New_Line;

        -- 速度
        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);

        i := i + 1;
    end loop;
end Ada0801;

xxxxxx@yyyyyy /Z
$ gnatmake Ada0801.adb

xxxxxx@yyyyyy /Z
$ Ada0801
0.01      0.00000      0.00000     0.48790    0.48790
0.02     48.45371     48.35571     0.97244    0.97146
0.03     48.12203     47.92670     1.45366    1.45073
0.04     47.79520     47.50319     1.93161    1.92576
0.05     47.47312     47.08509     2.40634    2.39661
0.06     47.15570     46.67226     2.87790    2.86333
0.07     46.84284     46.26460     3.34633    3.32598
0.08     46.53443     45.86201     3.81167    3.78460
0.09     46.23039     45.46436     4.27398    4.23924
省略
6.20      9.24481    -23.75680   125.68321    2.37768
6.21      9.22125    -23.79424   125.77543    2.13974
6.22      9.19771    -23.83152   125.86740    1.90143
6.23      9.17422    -23.86864   125.95915    1.66274
6.24      9.15076    -23.90561   126.05065    1.42368
6.25      9.12734    -23.94242   126.14193    1.18426
6.26      9.10395    -23.97907   126.23297    0.94447
6.27      9.08060    -24.01557   126.32377    0.70431
6.28      9.05728    -24.05191   126.41435    0.46379
6.29      9.03401    -24.08809   126.50469    0.22291
6.30      9.01076    -24.12412   126.59479   -0.01833

VB.NET

Option Explicit

Module VB0801
    '重力加速度
    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(1) As Double
        vx(0) = v * Math.Cos(radian)
        '鉛直方向の速度
        Dim vy(1) As Double
        vy(0) = v * Math.Sin(radian)
        '経過秒数
        Dim t As Double = 0.0
        '位置
        Dim x As Double = 0.0
        Dim y As Double = 0.0

        'Euler法
        Dim i As Integer = 1
        Do While (y >= 0.0)
            '経過秒数
            t = i * h

            '位置
            x += h * vx(0)
            y += h * vy(0)

            Console.WriteLine(String.Format("{0,4:F2}{5}{1,8:F5}{5}{2,9:F5}{5}{3,9:F5}{5}{4,8:F5}", t, vx(0), vy(0), x, y, vbTab))

            '速度
            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)

            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 VB0801.vb

Z:\>VB0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

C#

using System;

public class CS0801
{
    // 重力加速度
    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[2];
        vx[0] = v * Math.Cos(radian);
        // 鉛直方向の速度
        double[] vy = new double[2];
        vy[0] = v * Math.Sin(radian);
        // 経過秒数
        double t = 0.0;
        // 位置
        double x = 0.0;
        double y = 0.0;

        // Euler法
        for (int i = 1; y >= 0.0; i++)
        {
            // 経過秒数
            t = i * h;

            // 位置
            x += h * vx[0];
            y += h * vy[0];
            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, 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];
        }
    }

    // 空気抵抗による水平方向の減速分
    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 CS0801.cs

Z:\>CS0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

Java

import static java.lang.System.out;

public class Java0801 {

    // 重力加速度
    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[2];
        vx[0] = v * Math.cos(radian);
        // 鉛直方向の速度
        double[] vy = new double[2];
        vy[0] = v * Math.sin(radian);
        // 経過秒数
        double t = 0.0;
        // 位置
        double x = 0.0;
        double y = 0.0;

        // Euler法
        for (int i = 1; y >= 0.0; i++) {
            // 経過秒数
            t = i * h;

            // 位置
            x += h * vx[0];
            y += h * vy[0];
            out.println(String.format("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f", 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];
        }
    }

    // 空気抵抗による水平方向の減速分
    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 Java0801.java

Z:\>java Java0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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[2];
    vx[0] = v * cos(radian);
    // 鉛直方向の速度
    double vy[2];
    vy[0] = v * sin(radian);
    // 経過秒数
    double t = 0.0;
    // 位置
    double x = 0.0;
    double y = 0.0;

    // Euler法
    for (int i = 1; y >= 0.0; i++)
    {
        // 経過秒数
        t = i * h;
        cout << setw(4) << fixed << setprecision(2) << t << "\t";

        // 位置
        x += h * vx[0];
        y += h * vy[0];
        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    << "\t";
        cout << setw(8) << fixed << setprecision(5) <<  y    << endl;

        // 速度
        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];
    }
    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 CP0801.cpp
cp0801.cpp:

Z:\>CP0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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[2];
    vx[0] = v * cos(radian);
    // 鉛直方向の速度
    double vy[2];
    vy[0] = v * sin(radian);
    // 経過秒数
    double t = 0.0;
    // 位置
    double x = 0.0;
    double y = 0.0;

    // Euler法
    int i;
    for (i = 1; y >= 0.0; i++)
    {
        // 経過秒数
        t = i * h;

        // 位置
        x += h * vx[0];
        y += h * vy[0];

        printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", t, vx[0], vy[0], x, y);

        // 速度
        vx[1] = vx[0] + h * fx(vx[0], vy[0]);
        vy[1] = vy[0] + h * fy(vx[0], vy[0]);
        vx[0] = vx[1];
        vy[0] = vy[1];
    }
    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 OC0801 OC0801.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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[2];
    vx[0] = v * cos(radian);
    // 鉛直方向の速度
    double vy[2];
    vy[0] = v * sin(radian);
    // 経過秒数
    double t = 0.0;
    // 位置
    double x = 0.0;
    double y = 0.0;

    // Euler法
    for (int i = 1; y >= 0.0; i++)
    {
        // 経過秒数
        t = i * h;

        // 位置
        x += h * vx[0];
        y += h * vy[0];

        writefln("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f", 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];
    }
}

// 空気抵抗による水平方向の減速分
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 D0801.d

Z:\>D0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01832

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[2] float64
    vx[0] = v * math.Cos(radian)
    // 鉛直方向の速度
    var vy[2] float64
    vy[0] = v * math.Sin(radian)
    // 経過秒数
    var t float64 = 0.0
    // 位置
    var x float64 = 0.0
    var y float64 = 0.0

    // Euler法
    for i := 1; y >= 0.0; i++ {
        // 経過秒数
        t = float64(i) * h

        // 位置
        x += h * vx[0]
        y += h * vy[0]

        fmt.Printf("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n", t, vx[0], vy[0], x, y)

        // 速度
        vx[1] = vx[0] + h * fx(vx[0], vy[0])
        vy[1] = vy[0] + h * fy(vx[0], vy[0])
        vx[0] = vx[1]
        vy[0] = vy[1]
    }
}

// 空気抵抗による水平方向の減速分
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 GO0801.go

Z:\>8l -o GO0801.exe GO0801.8

Z:\>GO0801
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

Scala

object Scala0801 {

    // 重力加速度
    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

        // Euler法
        euler(1, vx, vy, x, y)
    }

    def euler(i:Int, vx:Double, vy:Double, x:Double, y:Double):Unit = {
        // 経過秒数
        val t = i * h

        // 位置
        val wx = x + h * vx
        val wy = y + h * vy
        println("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f".format(t, vx, vy, wx, wy))

        // 速度
        val wvx = vx + h * fx(vx, vy)
        val wvy = vy + h * fy(vx, vy)

        if (wy >= 0.0)
            euler(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 Scala0801.scala
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

F#

module Fs0801

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)

// Euler法
let rec euler(i:int) (vx:double) (vy:double) (x:double) (y:double):unit =
    // 経過秒数
    let t = float(i) * h

    // 位置
    let wx = x + h * vx
    let wy = y + h * vy

    printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t vx vy wx wy

    // 速度
    let wvx = vx + h * (fx vx vy)
    let wvy = vy + h * (fy vx vy)

    if wy >= 0.0 then
        (euler (i+1) wvx wvy wx wy)
    else
        ()

// Euler法
(euler 1 vx vy x y)

exit 0

Z:\>fsi  --nologo --quiet Fs0801.fs
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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))))

; Euler法
(defn euler[i vx vy x y]
    ; 経過秒数
    (def t (* i h))

    ; 位置
    (def wx (+ x (* h vx)))
    (def wy (+ y (* h vy)))

    (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t vx vy wx wy))

    ; 速度
    (def wvx (+ vx (* h (fx vx vy))))
    (def wvy (+ vy (* h (fy vx vy))))

    (if (>= wy 0.0)
        (euler (+ i 1) wvx wvy wx wy)
        nil))

(euler 1 vx vy x y)

Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0801.clj
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833

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)

-- Euler法
euler::Integer->Double->Double->Double->Double->IO ()
euler i vx vy x y =
    let
        t = (fromIntegral i) * h
        wx = x + h * vx
        wy = y + h * vy
        wvx = vx + h * (fx vx vy)
        wvy = vy + h * (fy vx vy)
    in
        if y >= 0.0
            then do
                printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t vx vy wx wy
                euler (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

    -- Euler法
    euler 1 vx vy x y

Z:\>runghc Hs0801.hs
0.01    48.79037     48.79037     0.48790    0.48790
0.02    48.45371     48.35571     0.97244    0.97146
0.03    48.12203     47.92670     1.45366    1.45073
0.04    47.79520     47.50319     1.93161    1.92576
0.05    47.47312     47.08509     2.40634    2.39661
0.06    47.15570     46.67226     2.87790    2.86333
0.07    46.84284     46.26460     3.34633    3.32598
0.08    46.53443     45.86201     3.81167    3.78460
0.09    46.23039     45.46436     4.27398    4.23924
省略
6.20     9.24481    -23.75680   125.68321    2.37768
6.21     9.22125    -23.79424   125.77543    2.13974
6.22     9.19771    -23.83152   125.86740    1.90143
6.23     9.17422    -23.86864   125.95915    1.66274
6.24     9.15076    -23.90561   126.05065    1.42368
6.25     9.12734    -23.94242   126.14193    1.18426
6.26     9.10395    -23.97907   126.23297    0.94447
6.27     9.08060    -24.01557   126.32377    0.70431
6.28     9.05728    -24.05191   126.41435    0.46379
6.29     9.03401    -24.08809   126.50469    0.22291
6.30     9.01076    -24.12412   126.59479   -0.01833