C++で常微分方程式(ルンゲ・クッタ法)

※再度修正
参考 http://www.enjoy.ne.jp/~k-ichikawa/Satellite3.html
初速 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[5]; 
    vx[0] = v * cos(radian); 
    // 鉛直方向の速度
    double vy[5]; 
    vy[0] = v * sin(radian); 
    // 経過秒数
    double t = 0.0;
    // 位置
    double x[5];
    x[0] = 0.0;
    double y[5];
    y[0] = 0.0;

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

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

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

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

        wx    = vx[0] + vx[3];
        wy    = vy[0] + vy[3];
        x[4]  = h *    wx;
        y[4]  = h *    wy;
        vx[4] = h * fx(wx, wy);
        vy[4] = h * fy(wx, wy);

        x[0]  += ( x[1] +  x[2] * 2 +  x[3] * 2 +  x[4]) / 6;
        y[0]  += ( y[1] +  y[2] * 2 +  y[3] * 2 +  y[4]) / 6;
        vx[0] += (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6;
        vy[0] += (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6;

        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.01	48.45619	 48.35852	  0.48623	 0.48574
0.02	48.12689	 47.93222	  0.96914	 0.96719
0.03	47.80236	 47.51133	  1.44878	 1.44440
0.04	47.48250	 47.09575	  1.92520	 1.91743
0.05	47.16722	 46.68537	  2.39845	 2.38633
0.06	46.85642	 46.28007	  2.86856	 2.85116
0.07	46.55001	 45.87974	  3.33559	 3.31195
0.08	46.24790	 45.48430	  3.79958	 3.76877
0.09	45.94999	 45.09363	  4.26056	 4.22165
0.10	45.65621	 44.70764	  4.71859	 4.67066
省略
6.20	 9.25063	-23.74826	125.74997	 2.41410
6.21	 9.22707	-23.78565	125.84236	 2.17643
6.22	 9.20355	-23.82289	125.93451	 1.93839
6.23	 9.18006	-23.85997	126.02643	 1.69997
6.24	 9.15661	-23.89690	126.11811	 1.46119
6.25	 9.13319	-23.93366	126.20956	 1.22203
6.26	 9.10982	-23.97027	126.30078	 0.98251
6.27	 9.08647	-24.00673	126.39176	 0.74263
6.28	 9.06317	-24.04303	126.48251	 0.50238
6.29	 9.03990	-24.07917	126.57302	 0.26177
6.30	 9.01667	-24.11516	126.66330	 0.02080
6.31	 8.99347	-24.15100	126.75335	-0.22053