数値的に解いた自由落下のプログラミング

Question

###前提・実現したいこと
時間ごとの位置を出力する自由落下のプログラムを作っているのですが、問題が二つあります。

①時間(t)=0の時、位置が0でないこと。
②0.01刻みで数値的に解いているはずなのですが、解析的に解いた数値と異なってしまっていること。

以上の二つです。

###該当のソースコード

#include <stdio.h>
#define g 9.8
#define v0 0
int main(void){
        double Y,y,t,dt,a,v,m;

        m=1.0;
        t=0.0;
        y=0.0;
        v=v0;
        a=0.0;

        dt=0.01;

        printf("t       y\n");

        for (t=0.0; t<2.0; t=t+dt){
                a = -g;
                v = v+a*dt;
                y = y+v*dt;
                Y=-g*t*t/2.0;
              
        printf("%f %f %f\n",t,y,Y);
        }
        return 0;
}

Comments

[kei344]

もしかして質問内容をエディタで作ってコピー＆ペースしていませんか？teratail のバグで記述内容が表示されないというものが以前有ったので、ペースト後に何文字か記入してみてください。

[ikedas]

t=0.0のとき位置が0でないのは、すでに漸化式を1回計算しているので当然だと思います (最初の結果はt = 0.01のときのものだということ)。なおできればでいいのですが、数値的に解いた数値の隣に解析的に解いた数値も並べて出力できませんか。何かわかるかもしれません。

Answer 1

1回目のループ(t=0.0のとき)にvとyが0.01秒のときを求めてしまっているからです。

v = v + a*dtは前回の速度から現在の速度を求める、y = y + v*dtは前回の位置と現在の速度から求めるですが、t=0.0のときを考えると、前回は-0.01秒前ですので、vもyも計算前の値は-0.01秒前の値でなければなりません。しかし、vとyの初期値はt=0.0の時のものであり、-0.01秒前はあらいません。t=0.0の時の値に対して、そのまま計算してしまうと次回の値、つまり、0.01秒後の値になってしまいますので、ずれが発生してしまいます。

解決方法としては、t=0.0の時のみ条件分岐して、前回の値から現在の値を求めるのでは無く、そのまま初期値を使えばうまく行くでしょう。なお、Yとは必ずずれていきますが、dtを小さくすると、誤差も小さくなっていくことが確認できると思います。