私はプログラム(ここではforranだが他の言語の方もアドバイスお願いします。)上での時間積分式にオイラー法を使用しているのですが精度が悪いのか途中で値がNaNになりうまくいきません。
そこでルンゲ・クッタ法など、より精度の良い時間積分を用いたいのですが調べてみてもプログラム上でどう実装するのかが分かりません。
今2次元の流体方程式を解こうとしていて、計算スキームは2次の中心差分法を用いています。↓のコードをオイラー法からより精度の良いものに変えたいです。長くなりましたが分かる方よろしくお願いします。
自分のコード(1部省略)# dnu/dt = R_nuの微分方程式を考える。C1= 1/(2×格子幅)
do ix = ~
do iy = ~
R_nu(ix,iy) = -ux(ix,iy)C1(nu(ix+1,iy)-nu(ix-1,iy))-uy(ix,iy)C1(nu(ix,iy+1)-nu(ix,iy-1)) &
-nu(ix,iy)C1(ux(ix+1,iy)-ux(ix-1,iy) +uy(ix,iy+1)-uy(ix,iy-1))
enddo
enddo
~省略~
nu(:,:) = nu(:,:) + R_nu*dt #時間発展式(オイラー法)
時間発展させると途中で値がNanになってしまいます
例）
#こんな感じでx座標 y座標 nu ux uy siitaというそれぞれ2次元変数に関する方程式を回して値を出しています。
-0.500000E+02 -0.500000E+02 0.121000E+00 NaN NaN 0.100000E+01
または、問題・エラーが起きている画像をここにドラッグアンドドロップ
例）
```ruby
def greet
puts Hello World
end
ここに問題・エラーに対して試したことを記載してください。
精度と安定性は別の問題です。
たしかにRKを使用することで多少「時間方向の」安定性は向上しますが、そもそも不安定の原因は移流方程式の勾配型
u_j \frac{\partial u_i}{\partial x_j}
をナイーブに離散化している点、つまり「空間方向の」差分が不安定なことだと思います。
ここで説明するには長くなりすぎるので結論だけ申し上げると、非圧縮性を仮定した上で発散型
\frac{\partial u_j u_i}{\partial x_j}
を離散化するのが最もわかりやすく簡単で効果があると思います。
現状だと解かれている式は流体方程式というよりは運動量に関する移流方程式ですが、非圧縮流体に対する主解法であるSMAC法やFractional step法について調べてみることをおすすめします。
Qiitaにも沢山記事がありますし、詳しく体系的に学びたい場合はこちらが良いと思います。
@hattorisena
Questioner