ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

C++で常微分方程式(中点法)

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

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

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

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

初期値 x0 から次の式によって, 順次 x1, x2, … を求める.

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

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

        // 位置・速度
        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;
        vx[0]     = vx[0] + h * fx(wx, wy);
        vy[0]     = vy[0] + h * fy(wx, wy);
        x[0]      =  x[0] + h *    wx;
        y[0]      =  y[0] + h *    wy;

        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.45620	 48.35854	  0.48622	 0.48573
0.02	48.12691	 47.93225	  0.96912	 0.96717
0.03	47.80240	 47.51138	  1.44876	 1.44438
0.04	47.48255	 47.09581	  1.92517	 1.91740
0.05	47.16729	 46.68543	  2.39841	 2.38629
0.06	46.85650	 46.28014	  2.86852	 2.85111
0.07	46.55009	 45.87983	  3.33554	 3.31189
0.08	46.24799	 45.48440	  3.79952	 3.76870
0.09	45.95010	 45.09374	  4.26050	 4.22158
0.10	45.65632	 44.70776	  4.71852	 4.67058
省略
6.20	 9.25071	-23.74815	125.75040	 2.41468
6.21	 9.22715	-23.78555	125.84279	 2.17702
6.22	 9.20362	-23.82279	125.93494	 1.93897
6.23	 9.18014	-23.85987	126.02686	 1.70056
6.24	 9.15669	-23.89679	126.11854	 1.46178
6.25	 9.13327	-23.93356	126.20999	 1.22262
6.26	 9.10989	-23.97017	126.30121	 0.98311
6.27	 9.08655	-24.00663	126.39219	 0.74322
6.28	 9.06325	-24.04293	126.48294	 0.50297
6.29	 9.03998	-24.07907	126.57346	 0.26236
6.30	 9.01674	-24.11506	126.66374	 0.02139
6.31	 8.99355	-24.15090	126.75379	-0.21994
参考文献