C++で常微分方程式(オイラー法)

初速 250 km/h で, 45°の角度で打ったボールの軌跡をオイラー法で計算する
(空気抵抗係数を 0.01 で計算)

重力による鉛直方向の減速分は, 重力加速度を g, 時間を t とすると,

空気抵抗による水平方向の減速分は,速度を v, 速度の水平方向成分を vx, 空気抵抗係数を k とすると,

同様に, 鉛直方向の減速分は, 速度の鉛直方向成分を vy とすると,

初期値 x₀ から次の式によって, 順次 x₁, x₂, … を求める.

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

// 重力加速度
const double g = -9.8;
// 空気抵抗係数
const double k = -0.01;
// 時間間隔(秒)
const double dt = 0.1;

// 空気抵抗による水平方向の減速分
double fx(double vx, double v);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v);

int main()
{
    // 角度 
    double degree = 45;
    double radian = degree * M_PI / 180.0;
    // 初速 250 km/h -> 秒速に変換
    double v = 250 * 1000 / 3600; 
    // 水平方向の速度
    double vx = v * cos(radian); 
    // 鉛直方向の速度
    double vy = 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 * dt;
        cout << setw(8) << fixed << setprecision(5) << t << "\t";

        // 速度
        cout << setw(8) << fixed << setprecision(5) << vx << "\t";
        cout << setw(9) << fixed << setprecision(5) << vy << "\t";

        // 位置
        x = x + dt * vx; 
        y = y + dt * vy; 
        cout << setw(9) << fixed << setprecision(5) << x << "\t";
        cout << setw(8) << fixed << setprecision(5) << y << endl;

        // 速度
        vx += fx(vx, v);
        vy += fy(vy, v);
        v  = sqrt(vx * vx + vy * vy);
    }
    return 0;
}

// 空気抵抗による水平方向の減速分
double fx(double vx, double v)
{
    double ax = k * v * vx;
    return ax * dt; 
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double vy, double v)
{
    double ay = k * v * vy;
    return (g + ay) * dt;
}

 0.10000	48.79037	 48.79037	  4.87904	 4.87904
 0.20000	45.42383	 44.44383	  9.42142	 9.32342
 0.30000	42.53716	 40.63944	 13.67514	13.38736
 0.40000	40.03469	 37.26862	 17.67861	17.11423
 0.50000	37.84493	 34.25015	 21.46310	20.53924
 0.60000	35.91324	 31.52194	 25.05442	23.69143
 0.70000	34.19713	 29.03567	 28.47413	26.59500
 0.80000	32.66301	 26.75310	 31.74044	29.27031
 0.90000	31.28395	 24.64356	 34.86883	31.73467
 1.00000	30.03808	 22.68214	 37.87264	34.00288
省略
 5.00000	12.23381	-18.02202	112.37476	27.65719
 5.10000	11.96733	-18.60946	113.57149	25.79624
 5.20000	11.70255	-19.17772	114.74174	23.87847
 5.30000	11.43964	-19.72687	115.88571	21.90578
 5.40000	11.17877	-20.25702	117.00358	19.88008
 5.50000	10.92013	-20.76834	118.09560	17.80325
 5.60000	10.66390	-21.26103	119.16199	15.67714
 5.70000	10.41025	-21.73532	120.20301	13.50361
 5.80000	10.15937	-22.19151	121.21895	11.28446
 5.90000	 9.91141	-22.62989	122.21009	 9.02147
 6.00000	 9.66655	-23.05081	123.17674	 6.71639
 6.10000	 9.42493	-23.45464	124.11924	 4.37093
 6.20000	 9.18669	-23.84177	125.03791	 1.98675
 6.30000	 8.95196	-24.21260	125.93310	-0.43451

参考文献