C++で常微分方程式(ホイン法)

初速 150 km/h で, 45°の角度で投げたボールの軌跡をホイン法で計算する
(空気抵抗を考慮しない)

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

// 重力加速度
const double g = 9.8;

// 水平方向の速度
double fx(double t, double vx)
{
    // 水平方向の速度は, 初速のまま   
    return vx; 
}
// 鉛直方向の速度
double fy(double t, double vy)
{
    // 初速から落下速度を引く 
    return vy - (t * g);
}

int main(void)
{
    // 角度 
    double degree = 45;
    double radian = degree * M_PI / 180.0;
    // 初速 150 km/h -> 秒速に変換
    double v = 150 * 1000 / 3600; 
    // 水平方向の速度
    double vx = v * cos(radian); 
    // 鉛直方向の速度
    double vy = v * sin(radian); 
    // 時間間隔(秒)
    double dt = 0.01;
    // 経過秒数
    double t = 0.0;
    // 位置
    double x = 0.0;
    double y = 0.0;

    // Heun法
    for (int i = 1; y >= 0.0; i++) 
    {
        // 位置
        // Euler法
//      x = x + dt * fx(t, vx);
//      y = y + dt * fy(t, vy);                             
       // Heun法
        x = x + (dt * fx(t, vx) + dt * fx(t + dt, vx)) / 2; 
        y = y + (dt * fy(t, vy) + dt * fy(t + dt, vy)) / 2; 
        // 経過秒数
        t = i * dt;

        // 経過秒数
        cout << setw(8) << fixed << setprecision(5) << t      << "\t";
        // 位置 (Heun法)
        cout << setw(8) << fixed << setprecision(5) << x      << "\t";
        cout << setw(8) << fixed << setprecision(5) << y      << "\t";
        // 位置 (正解)
        cout << setw(8) << fixed << setprecision(5) << vx * t << "\t";
        cout << setw(8) << fixed << setprecision(5) << vy * t - (g * t * t) / 2 << endl;
    }
    return 0;
}

 0.01000	 0.28991	 0.28942	 0.28991	 0.28942
 0.02000	 0.57983	 0.57787	 0.57983	 0.57787
 0.03000	 0.86974	 0.86533	 0.86974	 0.86533
 0.04000	 1.15966	 1.15182	 1.15966	 1.15182
 0.05000	 1.44957	 1.43732	 1.44957	 1.43732
 0.06000	 1.73948	 1.72184	 1.73948	 1.72184
 0.07000	 2.02940	 2.00539	 2.02940	 2.00539
 0.08000	 2.31931	 2.28795	 2.31931	 2.28795
 0.09000	 2.60922	 2.56953	 2.60922	 2.56953
 0.10000	 2.89914	 2.85014	 2.89914	 2.85014
省略
 5.80000	168.14999	 3.31399	168.14999	 3.31399
 5.81000	168.43991	 3.03502	168.43991	 3.03502
 5.82000	168.72982	 2.75506	168.72982	 2.75506
 5.83000	169.01973	 2.47412	169.01973	 2.47412
 5.84000	169.30965	 2.19221	169.30965	 2.19221
 5.85000	169.59956	 1.90931	169.59956	 1.90931
 5.86000	169.88948	 1.62544	169.88948	 1.62544
 5.87000	170.17939	 1.34058	170.17939	 1.34058
 5.88000	170.46930	 1.05474	170.46930	 1.05474
 5.89000	170.75922	 0.76793	170.75922	 0.76793
 5.90000	171.04913	 0.48013	171.04913	 0.48013
 5.91000	171.33904	 0.19135	171.33904	 0.19135
 5.92000	171.62896	-0.09840	171.62896	-0.09840

参考文献