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[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
参考文献