偏微分方程式: 差分法
偏微分方程式の数値解法: もう少し、踏み込む手法の一つ.有限差分法(finite difference method).
これまで何度か登場したが、「微分」を離散近似するシンプルな手法に、差分法と呼ばれるものがある. これは微分を、引き算と割り算で近似するもので、前回の method of line で右辺を離散近似する際でも使われていた.
この差分法を全面的に使って、微分方程式の微分をすべて差分で近似してしまって、偏微分方程式を完全に離散的な式(スキーム)で近似してしまおう、というのが今回紹介する方法だ.
この手法は単純に見えるが、単純なだけに応用は広いし、また、近似の仕方のバリエーションが豊富なため実は奥も深い. そしてなにより、うまくはまれば大変強力な数値解法になる.
今回、この手法についてまずは入門的な内容で様子をみてみよう.
対象は KdV 方程式で
前回紹介した KdV 方程式として, パラメータ $\epsilon$ をもつ次の式
\[ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0 \]
を考えることにしよう.
差分法を使ってみる
差分法は微分を差分、つまり、Taylor 展開で近似するものである. 近似できていれば良いのだ! と考えれば、こんなに自由な方法もないわけで、まあやりたい放題といって良い. 自由であるということと、良い近似解法ができるかどうかは全く別だが.
そこで、たとえば今回は、あちこち平均したり差分したりして作った以下の離散式(数値スキーム)、をもとの KdV 方程式の近似式として考えよう. ただし、$u_k^{(n)} \cong u(k\Delta x, n\Delta t)$ として、下の式は簡潔に書くために $u = u_k^{(n)}$, $u^{+} = u_k^{(n+1)}$ と記している.
\[ \frac{ u^{+} - u }{\Delta t} + \left( \frac{u^{+} + u}{2} \right) \delta_k^{(1)} \left( \frac{u^{+} + u}{2} \right) + \epsilon^2 \delta_k^{(3)} \left( \frac{u^{+} + u}{2} \right) = 0 \]
ただし、$n = 0,1,2, \cdots$, $k = 0, 1, \cdots, N-1$ で、離散境界条件は以前と同じように周期的境界条件を離散化した
\[ \left\{ \begin{array}{rcl} \displaystyle u_{-2} & := & u_{N-2}, \cr \displaystyle u_{-1} & := & u_{N-1}, \cr \displaystyle u_{N} & := & u_{0}, \cr \displaystyle u_{N+1} & := & u_{1} \cr \end{array} \right. \]
としよう.
早速計算してみよう!
では早速計算へ向けて準備を始めよう.
まずは wikipedia https://ja.wikipedia.org/wiki/KdV%E6%96%B9%E7%A8%8B%E5%BC%8F でも引用されている Zabusky and Kruskal の1965年の古い論文 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われているパラメータ等にあわせよう.
|
|
配列の添字を自由にしたいので、OffsetArrays package を用いるよ.
|
|
境界条件を適用するシーンが何回かありそうなので、関数にしておこう.
|
|
型を一回実際に使っておくと、あとで楽なので、初期値を作ってしまおう.まず、
|
|
としてから、Zabusky and Kruskal https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われている初期値を設定しよう.
|
|
OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}} with indices -2:201: 0.998027 0.999507 1.0
0.999507 0.998027 0.995562 0.992115 0.987688 0.982287 0.975917 0.968583 0.960294 0.951057 ⋮
0.951057 0.960294 0.968583 0.975917 0.982287 0.987688 0.992115 0.995562 0.998027 0.999507 1.0
0.999507
一階から三階まで、差分作用素を関数として定義しておこう.
|
|
数値スキームそのものを関数にしてしまおう
(すぐ後で書くように NLpackage で使う形式にするために) $0 = 式$ となるように移項して、この式部分を記述する格好にする. ただし、全体に $\Delta t$ をかけておくと良いかなあ. 式全体のオーダーを $u$ に合わせると、なにかと見易いのだ.
|
|
そして、この式がゼロになるような解 up を「探す」ことで時間ステップが進んだ値である up を求めるわけだが、それには以前も使った NLsolve package が使える(お手軽だね!).
しかし、NLsolve package は通常の配列しか知らないので、NLsolve package へそうした形で情報を渡せるように、通常の、外へはみ出てない配列と offsetarray を相互に変換する関数を NLsolve 用に作っておこう.
|
|
これを使って、数値スキームを NLsolve で使える形に直しておこう.
|
|
これで良い.
あとはこの f_solve 関数がゼロになるように数値解 up_in を探せ! と NLvole に命令すれば良いことになる. まずいつものように NLsolve を使う準備をしよう.
|
|
試しに 1ステップだけ解いてみよう.
|
|
200-element Array{Float64,1}: 0.999995 0.9996
0.998218 0.99585 0.992499 0.988167 0.982858 0.976578 0.969332 0.961127 0.951972 0.941874 0.930843 ⋮
0.928704 0.939882 0.950135 0.959453 0.967827 0.975247 0.981708 0.987201 0.991721 0.995264 0.997826 0.999404
図にしてみて、変でないことをたしかめておこう.
|
|
Plots.GRBackend()
|
|
ふむ、こんな感じか. 破綻等はしていないようで、まずは一安心.
注: Plots パッケージのバージョン等によっては、描かれる図がたいへん小さくなったりする.どうやら、画面の解像度の想定値がおかしい様子だ.そういうときは、
|
|
などとして、dpi を強制的に指定しておくと、それ以降はこの設定で描かれるのでそうしておくと良い.
さて、これで準備は整ったので、早速、計算してみよう
その前に
このスキームは非線形方程式を解くので、時間がかかることが想定される. そこで、計算がどれくらい進んでいるのか、リアルタイムで見えるようにしよう.それには、ProgressMeter というパッケージを使えば良い.
今の JuliaBox にはデフォルトでインストールされているので、JuliaBox 利用者はそのまま使える. 自分の PC などでは、 Pkg.add コマンドでインストールしておこう.
これを使うと、@showpprogress 時間間隔(秒) メッセージ などを for 文の頭につけるだけで、その for 文の実行状況と、残り推定時間を指定した時間間隔毎に更新して教えてくれる. 使い方はこれから見る例で十分だろう.
|
|
これでよし.では計算しよう. ただし、次の計算は juliabox の場合、サーバの混み具合によっては計算時間がたいへんかかる (通常の PC だと 1分程度なんだけどね). そうした場合は N や for 文の数を小さくするなど、現実的な計算時間になるようにパラメータを調整してやり直そう.
|
|
Progress: 100%|█████████████████████████████████████████| Time: 0:08:02
試しに結果の数字をちょろっと見てみる.
|
|
200×11 Array{Float64,2}: 1.0 0.956153 0.859704 0.757938 … -0.538929 -0.66376 -0.667301 0.999507 0.964282 0.871752 0.770684 -0.576519 -0.630957 -0.602904 0.998027 0.971663 0.88342 0.783294 -0.559205 -0.553665 -0.497812 0.995562 0.978272 0.89469 0.79573 -0.479927 -0.417903 -0.336871 0.992115 0.984084 0.905545 0.807932 -0.332199 -0.215148 -0.110155 0.987688 0.989074 0.915966 0.81983 … -0.100764 0.0718194 0.196364 0.982287 0.993214 0.925932 0.831366 0.221209 0.44255 0.57748 0.975917 0.996478 0.935423 0.842516 0.631325 0.881213 1.00676 0.968583 0.998838 0.944416 0.853315 1.08832 1.32207 1.41053 0.960294 1.00027 0.952887 0.86386 1.50778 1.66292 1.68633 0.951057 1.00073 0.960811 0.8743 … 1.76394 1.78613 1.73734 0.940881 1.00021 0.968161 0.884792 1.76168 1.64679 1.54699 0.929776 0.998671 0.974907 0.895442 1.49928 1.29714 1.18418 ⋮ ⋱ ⋮
0.929776 0.808463 0.690631 0.594284 0.412862 1.14712 1.32352 0.940881 0.823853 0.70618 0.608637 0.740933 1.32761 1.14066 0.951057 0.838747 0.721498 0.622884 … 1.04426 1.32333 0.835582 0.960294 0.853127 0.736578 0.63702 1.25063 1.13889 0.489555 0.968583 0.866973 0.751409 0.651037 1.29351 0.829242 0.160083 0.975917 0.880266 0.76598 0.664923 1.15866 0.477859 -0.11554 0.982287 0.892988 0.780279 0.678666 0.886113 0.143983 -0.33151 0.987688 0.905118 0.794297 0.692256 … 0.551663 -0.133731 -0.487255 0.992115 0.916636 0.80802 0.705686 0.220001 -0.349151 -0.59546 0.995562 0.92752 0.821436 0.718959 -0.0644842 -0.500781 -0.660936 0.998027 0.937749 0.834531 0.732083 -0.287025 -0.600866 -0.69433 0.999507 0.947301 0.847292 0.745072 -0.442895 -0.652234 -0.695196
あまり変なことにはなっていなさそうなので、プロットしてみよう.
|
|
|
|
|
|
以前同様、ソリトンへの分解が見られる.
ソリトンはその「高さ」によって移動速度が異なるため(背の高いほうが速い)、時間発展に従って速いものが前に素早く移動し、遅いものはだんだんとりのこされる.いま計算した結果は、だんだん分離するソリトンを見ていることになる.
前回も記したが、 wikipedia https://ja.wikipedia.org/wiki/KdV%E6%96%B9%E7%A8%8B%E5%BC%8F の右下に、同じ計算を motion gif にしたものが掲載されているので見てみよう.
今回も KdV 方程式の保存量をチェックしよう
前回も書いたように、KdV 方程式には、無限個の保存量があることがわかっていて、単純な方から挙げてみると、
\[ \left\{ \begin{array}{rcl} \displaystyle I_1 & = & \int_0^L u dx, \cr\cr \displaystyle I_2 & = & \int_0^L u^2 dx, \cr\cr \displaystyle I_3 & = & \int_0^L \left\{ u^3 / 3 - \epsilon^2 (u_x)^2 \right\}dx, \cr\cr \displaystyle & \vdots & \end{array} \right. \]
という感じだ. それぞれ、プログラムしてみよう.
|
|
早速、それぞれの量の時間発展を計算してみよう.
|
|
11-element Array{Float64,1}: -2.44249e-17 -2.05391e-16 -6.67244e-16 -1.56986e-15 -8.54872e-16 -7.64944e-16 -1.06914e-15 2.07057e-16 2.22045e-15 6.03517e-15 9.34253e-15
|
|
11-element Array{Float64,1}: 1.0
1.00003 1.00015 1.00082 1.00576 1.01573 1.02334 1.02699 1.02791 1.0273 1.02588
|
|
11-element Array{Float64,1}: -0.0047765
-0.00477688 -0.00477862 -0.00475651 -0.00395908 -0.002007
-0.000879984 -0.000662641 -0.000831951 -0.0011007
-0.00135988
とりあえず、この3つの保存量はだいたい保存されていることがみてとれる.
複数のソリトンの挙動を追いかけてみよう
先にも書いたように、ソリトンはピーク値が大きい(背が高い)ほうが速く移動する.
そこで、二つのソリトンを並べた初期値を用意して、速いほうが遅い方を追い越すシチュエーションを計算してみよう. このとき、phase shift とよばれる現象(追い越すほうが前へ、追い越される方が後ろへジャンプする)が発生するとされているので、それを観測してみよう.
(疑似) 2-ソリトン解
先にも示したように、高さ $H$, 中心 $x_0$ の 1-ソリトン解は次の式で書ける.
\[ u(x) = H sech^2( \sqrt{\frac{H}{12\epsilon^2}} (x-x_0) ) \]
そこでこれまた前回同様、少し離してこれを二つ配置することで 2-ソリトンを擬似的に初期配置してみよう.
|
|
OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}} with indices -2:201: 1.035e-6
8.59703e-7 8.00653e-6 1.04091e-5 1.35327e-5 1.75936e-5 2.28731e-5 2.97369e-5 3.86603e-5 5.02614e-5 6.53436e-5 8.49514e-5 0.000110443 ⋮
4.56736e-6 3.7938e-6
3.15126e-6 2.61754e-6 2.17421e-6 1.80597e-6 1.5001e-6
1.24603e-6 1.035e-6
8.59703e-7 8.00653e-6 1.04091e-5
試しにプロットして確認しておこう
|
|
では早速計算してみよう.
これまた、計算時間が膨大になりそうなときは、現実的な計算時間になるようにパラメータを調整してやり直そう.
(通常の PC だと 3分ぐらいなんだけどね…)
|
|
Progress: 100%|█████████████████████████████████████████| Time: 0:14:38
ではプロットしてみよう
|
|
ちょっとわかりにくいので上から見よう.計算パラメータが粗いので、少し図が「ガタついている」のはしかたないね.
|
|
下側の辺が初期値で、全体として上へ時間発展している図になる.
やはり、二つのソリトンが基本的に一定速度で動いていること、 背の高い、速いソリトンが後ろから追い越すときに、phase shift が起きていること、などがきちんと見て取れる.
この計算では保存量は?
|
|
26-element Array{Float64,1}: 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198
|
|
26-element Array{Float64,1}: 0.137543 0.137547 0.137548 0.137548 0.137547 0.137546 0.137543 0.137541 0.137537 0.13753 0.137513 0.137489 0.137449 0.137391 0.137318 0.137245 0.137205 0.137218 0.137277 0.137354 0.13742 0.137468 0.137502 0.13752 0.137532 0.137539
|
|
26-element Array{Float64,1}: 0.0239467 0.0239477 0.0239478 0.0239479 0.0239477 0.023947 0.023946 0.0239448 0.0239425 0.0239386 0.0239307 0.0239187 0.0238984 0.0238681 0.0238293 0.0237908 0.0237688 0.0237758 0.0238076 0.0238482 0.023883 0.0239082 0.023925 0.0239345 0.0239403 0.0239437
|
|
|
|
|
|
二つのソリトンが衝突するあたりで数値誤差が発生している様子など、先週の method of line + Runge-Kutta 法による結果とよく似ている.
Report No.10 差分法で計算
- 4つ目の保存量チェック
wikipedia の「保存量」項 を参照して4つ目の保存量 $I_4$ を調べて、これ(の近似値)が数値計算でどう変化するか、前回(method of line)や今回(差分法)の計算方法に対して調べてみよ.グラフも描こう. - 動画チャレンジ
上の方法で計算した KdV 方程式の数値解を動画にしよう. - 上の差分法スキームを、もっと高速に計算できないだろうか
いろいろ考えて試してみて、それぞれの計算時間を比較しよう.
ただし、計算結果があまり崩れるのも意味がなさそうなので、上の3つの保存量チェックぐらいはしておこう.