自由落下の方程式

まず、高校で、以下のような運動方程式を習います。 $$ ma = F $$ 自由落下では、質点が重力を受けることを考慮するので、重力加速度gを用いて、 $$ ma = -mg $$ 両辺から、mを除して、 $$ a = -g $$ 今回求めるのが、位置xだとして、 $$ \frac{d^{2}x}{{dt}^2}=-g $$ これが、解かなければならない微分方程式になります。
この方程式は、2階の情微分方程式なので、前に紹介したオイラー法で解くことができます。
オイラー法で解く際に、速度vを導入し、以下のように式を変形します。 $$ \frac{dx}{dt} = v $$ $$ \frac{dv}{dt} = -g $$ 上式を踏まえて、
オイラー法では、次のように、xとvを更新していきます。

$$ x_{i+1} = x_i + vdt $$

$$ v{i+1} = v{i} - gdt $$

自由落下を解くC言語コード

今回は、問題設定として、上向きに初期速度100で投げて、地面に落ちてくるまでをシミュレートします。
時間刻み$dt$、速度の初期値$v_0$、位置の初期値$x_0$を次のようにします。 $$ dt = 0.01 $$ $$ v_0 = 100.0 $$ $$ x_0 = 0.0 $$ これを踏まえて、厳密解は以下のように求まります。

$$ x(t) = {-0.5gt}^{2} + 100.0t $$

$$ v(t) = -gt + 100.0 $$ これを踏まえて、コードは以下のようになります。

#include <stdio.h>
#include <math.h>

double f_x(double x_pre, double v_pre, double dt);
double f_v(double x_pre, double v_pre, double g, double dt);
double calc_true_x_val(double t, double g);
double calc_true_v_val(double t, double g);

int main(void){
    FILE *outputfile;
    outputfile = fopen("out.txt", "w");
    // パラメータ
    double m = 1.0;
    double g = 9.80665;
    // 初期値
    double x = 0.0;
    double v = 100.0;
    // 解析範囲
    double t_start = 0.0;
    double t_end = 50.0;
    // 刻み時間
    double dt = 0.01;
    // 分割数
    int n = (int)((t_end-t_start) / dt) + 1;

    printf("t, x, v, true_x, true_v\n");
    for(double i=0; i<n; i++){
        double x_new = f_x(x, v, dt);
        double v_new = f_v(x, v, g, dt);
        double true_x_val = calc_true_x_val(t_start+dt*i, g);
        double true_v_val = calc_true_v_val(t_start+dt*i, g);
        x = x_new;
        v = v_new;
        printf("%5.5lf,%5.5lf,%5.5lf,%5.5lf,%5.5lf\n", t_start+dt*i, x, v, true_x_val, true_v_val);
        fprintf(outputfile, "%5.5lf,%5.5lf,%5.5lf,%5.5lf,%5.5lf\n", t_start+dt*i, x, v, true_x_val, true_v_val);
        if(x_new<=0.0){
            break;
        }
    }
}

double f_x(double x_pre, double v_pre, double dt){
    return x_pre + v_pre * dt;
}

double f_v(double x_pre, double v_pre, double g, double dt){
    return v_pre - g * dt;
}

double calc_true_x_val(double t, double g){
    return -0.5*g*t*t + 100.0*t;
}

double calc_true_v_val(double t, double g){
    return -g*t + 100.0;
}