緩募: Rust・Gnuplot・ffmpeg による 1 次元波動方程式の差分解法のバグ潰し¶
追記¶
次のツイートの指摘を受けて修正した.
差分化の仕方を変えれば挙動も当然変わってしまいます。 pic.twitter.com/EJLcgVdPP1
— 黒木玄 Gen Kuroki (@genkuroki) March 31, 2020
既に黒木氏からのコメントがありますが、まずvをアップデートした後にuをv(t,x)ではなくv(t+dt,x)でアップデートする必要があるようです。僕が間違っていなければ、いわゆるvon Neumannの安定性解析でv(t,x)だと解が必ず不安定になることが示せます（僕はみょんちゃんが言っていたので気づきました）
— Q (@life_wont_wait) March 31, 2020
これも結論を言うと, そもそもの離散化の仕方・計算の仕方がまずかった. 詳しくは上記黒木さんのツイートまたは GitHub 参照. Julia も気になっているが, まずは Rust でベクトル計算できないかが気になる. てらモスの兄貴に聞くのがいいだろうか?
結論から¶
結論から言うと, 次のコードが意図通りに動かない.
これが次のような挙動をする.
意図しているのは次のような挙動.
バグなしコードの説明¶
追記: 最初にも追記したように, 修正済. 詳しくは上記黒木さんのツイートまたは GitHub 参照.
まずはバグなしコードから.
このページを参考にして, 素直に波動方程式を解いている. 初期条件は $u = 0, u_t = 0$ で, 原点で波を $(1/2) \sin (2 \pi f t)$ で駆動していて, 次のように離散化している.
\begin{align} u(t + \Delta t, x) = 2 u(t,x) - u(t - \Delta t, x) + \left( \frac{\Delta t}{\Delta x} \right)^2 (u(t,x + \Delta x) - 2 u(t,x) + u(t,x)). \end{align}
意図通りに動いているので問題ないだろうと思っている.
バグありコードの説明¶
初期条件や波の駆動は 1dim_wave_eq_only_u.rs と同じ設定にしている. 2 階の常微分方程式でよくやるように, $u, v = u_t$ として, 次の連立の偏微分方程式を解く形にしている.
\begin{align} u_t = v, \quad v_t = u_{xx}. \end{align}
これを次のように離散化した.
[u(t + \Delta t, x) = u(t, x) + v(t, x) \Delta t,]
[v(t + \Delta t, x) = v(t, x) + \frac{u(t, x - \Delta x) - 2 u(t,x) + u(t, x + \Delta x)}{2 \Delta x} \Delta t.]
実際には次のように実装している.
どこがまずいのかがよくわかっていない. どなたかわかる方, ぜひ教えてほしい.