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

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

#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 dt, double vx, double v);
// 重力と空気抵抗による鉛直方向の減速分
double fy(double dt, 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(8) << fixed << setprecision(5) << vy << "\t";

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

        // 速度
        v = sqrt(vx * vx + vy * vy);
        float vx1 = vx + fx(dt, vx, v);
        float vy1 = vy + fy(dt, vy, v);
        for (int j = 0; j <= 10; j++)
        {
            v = sqrt(vx1 * vx1 + vy1 * vy1);
            float vx2 = vx + fx(dt, vx1, v);
            float vy2 = vy + fy(dt, vy1, v);
            if ((fabs(vx1 - vx2) < 0.00001) && (fabs(vy1 - vy2) < 0.00001)) break;
            vx1 = vx2;
            vy1 = vy2;
         }
        vx = vx1;
        vy = vy1;
    }
    return 0;
}

// 空気抵抗による水平方向の減速分
double fx(double dt, double vx, double v)
{
    double ax = k * v * vx;
    return ax * dt; 
}
// 重力と空気抵抗による鉛直方向の減速分
double fy(double dt, 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.84741	44.92653	  9.46378	 9.37169
 0.30000	43.25562	41.46219	 13.78934	13.51791
 0.40000	40.95800	38.33189	 17.88514	17.35110
 0.50000	38.90908	35.48337	 21.77605	20.89943
 0.60000	37.07220	32.87448	 25.48327	24.18688
 0.70000	35.41746	30.47085	 29.02502	27.23397
 0.80000	33.92023	28.24416	 32.41704	30.05838
 0.90000	32.56007	26.17089	 35.67305	32.67547
 1.00000	31.31983	24.23135	 38.80503	35.09861
省略
 5.00000	12.88588	-16.83596	116.88782	34.31772
 5.10000	12.61436	-17.44057	118.14926	32.57366
 5.20000	12.34465	-18.02671	119.38372	30.77099
 5.30000	12.07688	-18.59443	120.59141	28.91155
 5.40000	11.81120	-19.14380	121.77253	26.99717
 5.50000	11.54775	-19.67495	122.92731	25.02967
 5.60000	11.28670	-20.18802	124.05598	23.01087
 5.70000	11.02821	-20.68321	125.15880	20.94255
 5.80000	10.77242	-21.16075	126.23604	18.82648
 5.90000	10.51948	-21.62089	127.28799	16.66439
 6.00000	10.26955	-22.06392	128.31494	14.45800
 6.10000	10.02277	-22.49016	129.31722	12.20898
 6.20000	 9.77926	-22.89993	130.29514	 9.91899
 6.30000	 9.53915	-23.29361	131.24906	 7.58963
 6.40000	 9.30255	-23.67155	132.17931	 5.22247
 6.50000	 9.06956	-24.03415	133.08627	 2.81906
 6.60000	 8.84029	-24.38180	133.97030	 0.38088
 6.70000	 8.61481	-24.71493	134.83178	-2.09062

参考文献