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 dt = 0.1; // 空気抵抗による水平方向の減速分 double fx(double vx, double v); // 重力と空気抵抗による鉛直方向の減速分 double fy(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(9) << fixed << setprecision(5) << vy << "\t"; // 位置 x = x + dt * vx; y = y + dt * vy; cout << setw(9) << fixed << setprecision(5) << x << "\t"; cout << setw(8) << fixed << setprecision(5) << y << endl; // 速度 vx += fx(vx, v); vy += fy(vy, v); v = sqrt(vx * vx + vy * vy); } return 0; } // 空気抵抗による水平方向の減速分 double fx(double vx, double v) { double ax = k * v * vx; return ax * dt; } // 重力と空気抵抗による鉛直方向の減速分 double fy(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.42383 44.44383 9.42142 9.32342 0.30000 42.53716 40.63944 13.67514 13.38736 0.40000 40.03469 37.26862 17.67861 17.11423 0.50000 37.84493 34.25015 21.46310 20.53924 0.60000 35.91324 31.52194 25.05442 23.69143 0.70000 34.19713 29.03567 28.47413 26.59500 0.80000 32.66301 26.75310 31.74044 29.27031 0.90000 31.28395 24.64356 34.86883 31.73467 1.00000 30.03808 22.68214 37.87264 34.00288 省略 5.00000 12.23381 -18.02202 112.37476 27.65719 5.10000 11.96733 -18.60946 113.57149 25.79624 5.20000 11.70255 -19.17772 114.74174 23.87847 5.30000 11.43964 -19.72687 115.88571 21.90578 5.40000 11.17877 -20.25702 117.00358 19.88008 5.50000 10.92013 -20.76834 118.09560 17.80325 5.60000 10.66390 -21.26103 119.16199 15.67714 5.70000 10.41025 -21.73532 120.20301 13.50361 5.80000 10.15937 -22.19151 121.21895 11.28446 5.90000 9.91141 -22.62989 122.21009 9.02147 6.00000 9.66655 -23.05081 123.17674 6.71639 6.10000 9.42493 -23.45464 124.11924 4.37093 6.20000 9.18669 -23.84177 125.03791 1.98675 6.30000 8.95196 -24.21260 125.93310 -0.43451
参考文献