計算結果の検証

　得られた計算結果は正しいのだろうか？

前ページのグラフをみると，よく知っている放物線軌道を描いており，良さそうな気がする． このように正解を知っているのならば，それと比べれば計算結果の検証ができる． しかし，そもそも正解を知っている問題をわざわざ数値計算で解く必要はない． 実際に数値計算が扱うのは手計算では解が得られない問題が大半である． このような場合には，どのようにして計算結果が正しいと判断すべきだろうか．

解析解との比較

ボール投げ問題に関しては解析解が分かっているので，まずはそれと比較してみよう． ボール投げの運動方程式の解析解は，

(1)\[\begin{split}x &= x_0 + v_0\cos{\theta}\ t \\ y &= y_0 + v_0\sin{\theta}\ t - \frac{1}{2}gt^2\end{split}\]

であった．先に示したプログラム >> の34行目から下を以下のように書き換えよう．

xa= x(1) + u(1)*t
ya = y(1) + v(1)*t - 0.5d0*g*t**2

open(unit=10,file='output.dat',status='replace')
write(10,*)'t, x, y, xa, ya'
do i = 1,n
    write(10,*)t(i),',',x(i),',',y(i),',',xa(i),',',ya(i)
end do
close(10)

Warning

この変更を素直に加えた後に，コンパイルを行ってエラーがでた人は．．．新たに使ってる配列xa,yaの宣言（double precision）を忘れずに．

xa= x(1) + u(1)*t
ya = y(1) + v(1)*t - 0.5d0*g*t**2

この2行で，新たに解析解の計算を行っている． xa,ya,tは配列であるが，このように配列全体に対して一括で演算することができる． 配列tの各要素に対して，u(1)とtを掛け算した結果にx(1)を足し算して，配列xaに格納している． 当然xa,yaとtの配列の大きさ（要素の数）は一致していなくてはいけない． もちろん，オイラー法の計算部分のように1番目の配列から順番に逐一計算を行っていってもよいが， 中身が全てわかっている配列（ここではt）を使って計算する場合は，この書き方のほうが楽である．

open(unit=10,file='output.dat',status='replace')
write(10,*)'t, x, y, xa, ya'
do i = 1,n
    write(10,*)t(i),',',x(i),',',y(i),',',xa(i),',',ya(i)
end do
close(10)

これまでは，計算結果をターミナルに表示していたが，情報も増えてきたのではじめからファイルに出力することにした．

ファイルに結果を出力する場合，まずはopen文によって，出力先のファイルを指定し，開く．

unit=10は，出力先の装置番号を意味する．何番でもよいが，5番はキーボード，6番はディスプレイと決まっていることが多いので，5,6以外の番号がよい．

file=’output.dat’はファイル名である．

status=’replace’は指定されたファイルが存在していない場合に新しくファイルを生成する．存在している場合には，新しファイルで置き換える．

status=’old’は指定されたファイルが存在していなくてはいけない．既存のファイルからデータを読み込む場合に主に使う．

status=’new’は指定されたファイルが存在していてはいけない．存在しているとエラーになる．

これまでwrite(＊,＊)としていた部分がwrite(10,*)に代わっている． この10はopenで生成したファイルの装置番号である．元々の*は標準出力を意味していてディスプレイ（6）になる．

Fig. Euler implicit and analytical solution

新たに計算した解析解を書き加えた．オイラー法による計算結果は確かに解析解と同じく放物線を描いており， その数値もほぼ一致している．しかし，完全に一致しているわけではなく差がある．またこの差は時間の経過とともに徐々に大きくなっているのが分かる．

保存量での検証

計算対象に保存量があるならば，それを計算するのがよいだろう． 保存量とは時間によって変化しない，一定を保つ量のことである． エネルギーを散逸させるような，摩擦とか粘性がない系ならば，全エネルギーが保存しているはずである． つまり，ボール投げ問題では，運動エネルギーとポテンシャル（位置）エネルギーの和が一定である．

(2)\[E = \frac{1}{2} m v^2 + mgh = constant\]

さきほど変更を加えた部分をさらに以下のように変更して全エネルギーの計算を行う．

xa= x(1) + u(1)*t
ya = y(1) + v(1)*t - 0.5d0*g*t**2
ua = u(1)
va = v(1) - g*t
E = 0.5d0*(u**2 + v**2) + g*y
Ea = 0.5d0*(ua**2 + va**2) + g*ya

open(unit=10,file='outputE.dat',status='replace')
write(10,*)'t, x, y, xa, ya, E, Ea'
do i = 1,n
    write(10,*)t(i),',',x(i),',',y(i),',',xa(i),',',ya(i),',',E(i),',',Ea(i)
end do
close(10)

質量mも一定なのでE/mを計算している． 解析解による運動エネルギー（Ea）は一定だと分かっているので，いちいち計算する必要はないが， 一応計算した．一定値である．

Fig. Total energy

オイラー法による結果では，全エネルギーが時間の経過とともに徐々に大きくなってしまっている．

ボールの軌道を比較したときには，解析解，数値解ともに放物線だったしまあいいかという気になった人もいるかもしれない． あるいは保存するはずの全エネルギーが増えてしまうのはおかしいと思った人もいるかもしれない．

Note

保存量がたしかに保存しているかどうかは，数値解析の精度の検証手段として有効である．ただし，保存量が保存しているからと言って，正しい結果が得られているとは限らない．