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

    // Heun法
    for (int i = 1; y[0] >= 0.0; i++) 
    {
        // 経過秒数
        t = i * h;
        cout << setw(5) << fixed << setprecision(3) << t << "\t";

        // 位置・速度
        x[1]  =  x[0] + h *    vx[0];
        y[1]  =  y[0] + h *    vy[0];
        vx[1] = vx[0] + h * fx(vx[0], vy[0]);
        vy[1] = vy[0] + h * fy(vx[0], vy[0]);

        x[2]  =  x[0] + h * (  vx[0]          +    vx[1]        ) / 2;
        y[2]  =  y[0] + h * (  vy[0]          +    vy[1]        ) / 2;
        vx[2] = vx[0] + h * (fx(vx[0], vy[0]) + fx(vx[1], vy[1])) / 2;
        vy[2] = vy[0] + h * (fy(vx[0], vy[0]) + fy(vx[1], vy[1])) / 2;

        x[0]  =  x[2];
        y[0]  =  y[2];
        vx[0] = vx[2];
        vy[0] = vy[2];

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

0.010	48.45620	 48.35853	  0.48622	 0.48573
0.020	48.12690	 47.93224	  0.96912	 0.96717
0.030	47.80238	 47.51136	  1.44876	 1.44437
0.040	47.48254	 47.09579	  1.92517	 1.91740
0.050	47.16726	 46.68541	  2.39841	 2.38629
0.060	46.85647	 46.28011	  2.86852	 2.85111
0.070	46.55007	 45.87980	  3.33554	 3.31189
0.080	46.24796	 45.48436	  3.79952	 3.76870
0.090	45.95006	 45.09369	  4.26050	 4.22158
0.100	45.65629	 44.70771	  4.71852	 4.67057
省略
6.200	 9.25067	-23.74819	125.75002	 2.41403
6.210	 9.22711	-23.78559	125.84241	 2.17636
6.220	 9.20359	-23.82283	125.93456	 1.93832
6.230	 9.18010	-23.85991	126.02648	 1.69991
6.240	 9.15665	-23.89683	126.11816	 1.46112
6.250	 9.13323	-23.93360	126.20961	 1.22197
6.260	 9.10985	-23.97021	126.30083	 0.98245
6.270	 9.08651	-24.00667	126.39181	 0.74256
6.280	 9.06321	-24.04297	126.48256	 0.50232
6.290	 9.03994	-24.07911	126.57307	 0.26171
6.300	 9.01671	-24.11510	126.66335	 0.02073
6.310	 8.99351	-24.15094	126.75341	-0.22060

参考文献