しぇんゆバッテリー - えくせんこりしてぃ

#ラム乾ありがとうございました！今週は謎を作れませんでした, すみません. 代わりと言っちゃなんですが第4話「ふたりの射程」をテーマに物理()の問題を作りました！

雨嘉も神琳もよかったね......#アサルトリリィ #アサルトリリィ夢結 #アサルトリリィ弾丸 #アサルトリリィ木曜日 pic.twitter.com/COv7Zyc4NC
— いのこり (@ImCallinx2You) 2020年10月29日

ってことで

ゆーじあさんの弾丸の軌道(第4話)を数値計算しました

誤差5cmで神琳撃てた
ここまで精度良ければ私が雨嘉でしょ pic.twitter.com/7cnhkFvBQp
— いのこり (@ImCallinx2You) 2020年11月13日

なわけないでしょ

神琳から雨嘉方向を見たときの弾道をプロットしました（時間が進むにつれて赤くなります）
パラメータに依るんですけど風で浮いてからギューンと落ちてきます#アサルトリリィ pic.twitter.com/E98hqUIwoI
— いのこり (@ImCallinx2You) 2020年11月27日

きれいだよ神琳

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0; i<(n); i++)
using namespace std;
using vd = vector<double>;
using vv = vector<vd>;

const double l  = 1000.0;
const double v0 = 1800.0;
const double g  =    9.8;
const double k  = 1e-3;
const double dt = 1e-3;
const double w0 = -10;

#define XLEN 6

vv init_f(double theta, double phi){
    vv x0(2, vd (XLEN, 0.0));
    
    x0[0][0] = 0;
    x0[0][1] = 0;
    x0[0][2] = 0;
    x0[0][3] = v0*cos(theta);
    x0[0][4] = v0*sin(theta)*sin(phi);
    x0[0][5] = v0*sin(theta)*cos(phi);
    
    return x0;
}

vd f(vd bx){
    vd dxdt(XLEN, 0);
    
    double v2 = pow(bx[3], 2) + pow(bx[4]-w0, 2) + pow(bx[5], 2);
    
    dxdt[0] = bx[3];
    dxdt[1] = bx[4];
    dxdt[2] = bx[5];
    dxdt[3] = -k*sqrt(v2)*bx[3];
    dxdt[4] = -k*sqrt(v2)*(bx[4]-w0);
    dxdt[5] = -k*sqrt(v2)*bx[5] - g;
    
    return dxdt;
}

vd rk4(vd bx, double dt){
    vd ax = bx, k, a = bx, h = {1,2,2,1,1};
    rep(j,4){
        k = f(a);
        rep(i,XLEN){
            k[i] *= dt;
            a[i] = bx[i] + k[i]/h[j+1];
            ax[i] += k[i]*h[j]/6.0;
        }
    }
    return ax;
}

vd shot(double theta, double phi, int n=0){
    vd yz(2, 0);
    vv x = init_f(theta, phi);
    double t = 0.0;
    while(x[0][0]<l){
        t += dt;
        x[1] = rk4(x[0], dt);
        if(n){
            printf("%lf ", t);
            rep(i,XLEN)printf("%.10lf%c", x[1][i], i==XLEN-1?'\n':' ');
        }
        swap(x[0],x[1]);
    }
    
    double a = x[0][0]-l;
    double b = l-x[0][1];
    
    rep(i,2)yz[i] = (x[1][i+1]*a+x[0][i+1]*b)/(a+b);
    
    return yz;
}

int main(){
    double vy_max = 0.1;
    double vy_min = 0.0;
    double vz_max = 0.1;
    double vz_min = 0.0;
    
    double r;
    double vy, vz;
    double theta, phi;
    vd mark;
    
    while(r>1e-3){
        
        vy = (vy_max+vy_min)/2.0;
        vz = (vz_max+vz_min)/2.0;
        theta = atan(sqrt(vy*vy + vz*vz));
        phi   = atan2(vy, vz);
        mark  = shot(theta, phi);
        
        if(mark[0]>0)vy_max = vy;
        else vy_min = vy;
        if(mark[1]>0)vz_max = vz;
        else vz_min = vz;
        
        r = sqrt(pow(mark[0], 2) + pow(mark[1], 2));
    }
    printf("\ntheta=%lf[deg], phi=%lf[deg]\n", theta*180/M_PI, phi*180/M_PI);
    printf("y=%lf[m], z=%lf[m]\n\n", vy*l, vz*l);
    
    return 0;
}

解くべき運動方程式
空気抵抗と風、重力が働くときの弾丸の軌道を計算します。
弾丸のように速度が大きい時、空気抵抗の大きさは二乗に比例します。これは弾丸前後面での空気の圧縮・希薄による圧力差が原因になっています*1。風がある時は、向かい風や追い風（=相対風速 ${\bf v}-{\bf w}_0$ ）をもとに空気抵抗が決まります。
したがって運動方程式は

$\displaystyle{\frac{d {\bf v}}{dt} = - k |{\bf v}-{\bf w}_0|\left({\bf v}-{\bf w}_0\right) + {\bf g}}$

で表されます。
ただし、毎秒1800 m という弾丸の速度は音速を超えているので二乗則が正常に働いているかどうかは分かりません......（文献を当たりましたが、雨嘉さんに撃てるほど小さな質量でのデータが見当たりませんでした）。

運動方程式の右辺第1項が風場ありの空気抵抗、第2項が重力を表します。空気抵抗と重力があるときの軌道は特殊な例を除いて解析解を持ちません。本問もそうなので、数値的に解を求めましょう。

コードの概説
まずアルゴリズムのフローチャートはこんな感じです。
試し撃ちを行い、着弾位置から判断して発射角度を修正することを繰り返し行います。雨嘉さんは一発で当てちゃうのですごい
f:id:Callinx2You:20210114184958j:plain

ではコードの解説に移ります。

個人的な読みやすさのために以下4つを使います。
rep(i,n): iを0からn-1まで行うforループ
vd: 実数の1次元配列
vv: 実数の2次元配列
XLEN: 6 (= 3次元 * 2)

コードは以下5つの関数から構成されます。

init_f：初期値の設定
f：解くべき運動方程式の指定
rk4：4次のルンゲ・クッタ法を用いた数値積分コア
shot：弾丸の軌道計算
main：神琳さんに命中するような角度を二分探索する

まず

init_f：初期値の設定

位置と速度の初期値 x0 を与えます。0~2が位置、3~5が速度に対応します*2。雨嘉さんの位置から $(\theta, \phi)$ 方向に速さ v0 で発射します。

vv init_f(double theta, double phi){
    vv x0(2, vd (XLEN, 0.0));
    
    x0[0][0] = 0;      // x
    x0[0][1] = 0;      // y
    x0[0][2] = 0;      // z
    x0[0][3] = v0*cos(theta);               // v_x
    x0[0][4] = v0*sin(theta)*sin(phi);      // v_y
    x0[0][5] = v0*sin(theta)*cos(phi);      // v_z
    
    return x0;
}

$\theta$ はx軸と速度方向のなす角、 $\phi$ はy軸方向への回転角を表します.

f：解くべき運動方程式の指定

運動方程式の右辺を記述します。

vd f(vd bx){
    vd dxdt(XLEN, 0);
    
    double v2 = pow(bx[3], 2) + pow(bx[4]-w0, 2) + pow(bx[5], 2);
    
    dxdt[0] = bx[3];      // dx/dt = v_x
    dxdt[1] = bx[4];      // dy/dt = v_y
    dxdt[2] = bx[5];      // dz/dt = v_z
    dxdt[3] = -k*sqrt(v2)*bx[3];           // dv_x/dt = a_x
    dxdt[4] = -k*sqrt(v2)*(bx[4]-w0);      // dv_y/dt = a_y
    dxdt[5] = -k*sqrt(v2)*bx[5] - g;       // dv_z/dt = a_z
    
    return dxdt;
}

状態 bxを使って微分係数 dxdtを返します。

rk4：4次のルンゲ・クッタ法を用いた数値積分コア

現在の状態 bxを使って、時間 dt後の状態 axを4次精度で計算します。

vd rk4(vd bx, double dt){
    vd ax = bx, k, a = bx, h = {1,2,2,1,1};
    rep(j,4){
        k = f(a);      // k_1~k_4を方程式fに基づいて順次計算
        rep(i,XLEN){
            k[i] *= dt;      // 刻み幅をかける
            a[i] = bx[i] + k[i]/h[j+1];      // 刻み幅の係数をかけて次の引数とする
            ax[i] += k[i]*h[j]/6.0;      // 戻り値に加算
        }
    }
    return ax;
}

ルンゲ＝クッタ法 - Wikipedia
めちゃめちゃ短く書けたから見て欲しい。
~~あのルンゲクッタさんをこんなに可愛くしちゃうなんてすごいことよ~~(これ5話だった)

shot：弾丸の軌道計算

$t=0$ から始めて、弾丸が $x=l$ を通過するまで数値積分を繰り返します。そして $x=l$ での $y$ , $z$ の値を通過前後の2点から線形補間して返します。

vd shot(double theta, double phi, int n=0){
    vd yz(2, 0);
    vv x = init_f(theta, phi);      // 与えられたtheta, phiで初期値を呼ぶ
    double t = 0.0;
    while(x[0][0]<l){               // xがlを超えたら終了
        t += dt;
        x[1] = rk4(x[0], dt);       // 次の時刻での位置・速度をルンゲクッタで計算
        if(n){      // (軌道出力用オプション)
            printf("%lf ", t);
            rep(i,XLEN)printf("%.10lf%c", x[1][i], i==XLEN-1?'\n':' ');
        }
        swap(x[0],x[1]);            // 状態の更新
    }
    
    double a = x[0][0]-l;
    double b = l-x[0][1];
    
    rep(i,2)yz[i] = (x[1][i+1]*a+x[0][i+1]*b)/(a+b);      // x=lでの通過位置を求める
    
    return yz;
}

main：神琳さんに（誤差 1 mmで）命中するような初期値を探索する

風と重力の影響は雨嘉さんからみて弾丸を右下に流すよう働くので、いい感じに左上に向けて撃たねば命中しません。そこで射撃シミュレーションを繰り返し行って、初速のy, z成分 ( $v_{0y}$ , $v_{0z}$ ) を二分探索します。
$x=l$ での弾丸の通過するy, z座標 ( $y_l$ , $z_l$ ) それぞれについて正負を判定します。
YES (通過位置が左/上) ならもっと初速を右/下にしてよいです。
NO (通過位置が右/下) なら初速は今より左/上でなければいけないです。

このアルゴリズムは、日本語で上手く言えないですが $\displaystyle{\frac{\partial y_l}{\partial v_{0z}}\sim 0}$ , $\displaystyle{\frac{\partial z_l}{\partial v_{0y}}\sim 0}$ を暗に仮定しています。
つまり通過位置が左だったからと言って最適初速がより右にあるとは限らないのです。
今回はたまたま上手くいったので良いですが、二分探索でこぼしてしまい解に収束しない可能性もある判定法でした。代わりに最急降下法などを使う方が自然です。

神琳さんから誤差 1 mm 以内になれば終了します。1周目での誤差が 100 mくらいあるので、20回程度のシミュレーションで収束することが期待されます。

int main(){
    
    //初速方向の探索範囲
    double vy_max = 0.1;
    double vy_min = 0.0;
    double vz_max = 0.1;
    double vz_min = 0.0;
    
    double r;
    double vy, vz;
    double theta, phi;
    vd mark;
    
    while(r>1e-3){            // 1 mm 以内になるまで繰り返す
        
        vy = (vy_max+vy_min)/2.0;
        vz = (vz_max+vz_min)/2.0;
        
        // vy, vz から theta, phi に変換
        theta = atan(sqrt(vy*vy + vz*vz));
        phi   = atan2(vy, vz);
        
        // 射撃シミュレーションを行う
        mark  = shot(theta, phi);
        
        // 探索範囲の更新
        if(mark[0]>0)vy_max = vy;
        else vy_min = vy;
        if(mark[1]>0)vz_max = vz;
        else vz_min = vz;
        
        // 神琳さんとの距離を求める
        r = sqrt(pow(mark[0], 2) + pow(mark[1], 2));
    }
    printf("\ntheta=%lf[deg], phi=%lf[deg]\n", theta*180/M_PI, phi*180/M_PI);
    printf("y=%lf[m], z=%lf[m]\n\n", vy*l, vz*l);
    
    return 0;
}

おわり。パラメータの説明や結果・考察はコピ本に載せました。

これは神雨だ......誰が何を言おうと神雨なんだ......

*1: 抗力 - Wikipedia 圧力抗力の項. 一方速度が小さい時は空気との摩擦が支配的で, 大きさは一乗に比例します(摩擦抗力).

*2:2行目は数値積分結果の更新のためにとってあります.

*3:そういえばrk4はt依存に対応していませんが, こちらのアップデートは予定していません. 楽そうだからやってもいいけど.