計算数理A演習第11回 – Masashi Shiraishi Homepage

１次元拡散方程式の数値解法（その２）

次の問題を考えます．先週の問題１．１と少し違いますので注意してください（初期値と x の範囲が異なる）．

問題　１．２

(1) $\begin{equation*} \left\{ \begin{array}{lc} \frac{\partial u(x, t)}{\partial t} = \frac{\partial^2u(x, t)}{\partial x^2} & (x \in (0, \pi), t> 0) \\ u(0, t) = 0 & (t>0) \\ u(\pi, t) = 0 & (t>0) \\ u(x, 0) = \sin(x) + \sin(7x) & ( x \in [0, \pi] ) \end{array} \right. \end{equation*}$

課題２

問題１．２の解を数値的に求めるプログラムを作成してください．

問題１．２の数値計算結果例

課題３

（初期値のsinの部分をcosに書き換えて）境界条件をノイマン条件に変更すると、どのような振る舞いをするか調べよ。
ノイマン条件にするには、式(4)の代わりに

$U_0^k = U_1^k, U_{N+1}^k = U_N^k \ \ (k=0, 1, 2, \cdots)$

とすればよい。

先に出てきた，

$U_j^{k+1} = \alpha U_{j-1}^k + (1-2\alpha)U_j^k + \alpha U_{j+1}^k$

を繰り返し計算し，近似解を求める方法を陽解法と呼びます．
（前の時刻 $t_n$ での値から次の時刻 $t_{n+1}$ の値を直ちに計算できる形の差分方程式を陽公式(explicit scheme)と呼びます．陽公式を使った解法なので陽解法と呼びます．）

天下り的になりますが，拡散方程式の数値解法として陽公式を用いる場合， $\tau$ と $h$ と $D$ は次の関係を満たす必要があります（安定性条件）．

(S) $\begin{equation*} \alpha = D\tau/h^2 < 1/2 \end{equation*}$

ここでは，この安定性条件を破った場合に計算がどのようになるかを見ます．

課題４

課題２のプログラムにおいてのパラメータ n （空間刻み数）やパラメータ　Ｍ（時間刻み数）　等を調節し，安定性条件（Ｓ）をわずかに破って見よ．

課題４に従い，安定性条件（Ｓ）をわずかに破った場合，次のような結果が得られます．
（Ｓ）の左辺の値が 0.5066 と，1/2 をわずかに越えているにすぎませんが，最終的には真の解からはほど遠いギザギザの状態になってしまいます．
1/2 を大きく越えると，計算がおかしくなるのも早まります．
（ここでは n および　M をα>1/2 となるように調整した．）

上の，ギザギザの解は，直感的に考えても熱方程式の解としてはおかしいことがわかると思います．（スプーンの熱分布がこの様になるとは思えません．）
陽解法においては安定性条件を満たさないパラメータで計算すると，計算が不安定となり，上記のような現象が見られます．
（ただし，これは，差分方程式の解としては正しいものです．しかし，近似したい熱方程式の近似解とはなっていところが問題なのです．）

注意：上の計算結果では対称性の崩れ（初期値は π/2 で対称）がみられます．これは，πを近似値として与えたためであると思われます．

この様に，陽解法においては，計算が安定に進むためには時間刻み幅および空間刻み幅が安定性条件(S)を満たさねばなりません．
ところが，この条件は計算量の点で問題があります．
例えば，時刻０から１まで計算を行うとします．
時間方向の格子点の総数は 1/τとなります．
計算の手間（計算時間）を減らす為には τ をなるべく大きく取りたい（τ が大きければ少ない繰り返し計算で目的の時刻まで計算できるので）．
そこで，安定性条件(S)より，

τ ≒ h²/(2D)

とすれば良いでしょう（(S)を余裕をもって満たす τ を用いる）．
一方，D は一定とした場合，精度の良い計算をするために空間に関する分解能を高めたいとします．
そこで h を 1/10 倍にしたとします．すると，τ は上の関係から 1/100 倍せざるを得ません．
ということは，時間方向の格子点数が 100倍に増えることになります．
これでは， h を小さくしていったときに，あまりに格子点数が多くなりすぎて計算の効率が悪くなります．
例えば今の例では，h を 1/10 倍したため，計算量は 10 倍となり（これは仕方が無い），更に安定性条件(S)を満たすために τ を 1/100 倍したので計算量は更に100倍となりました．
つまり，もとの計算量に比べて 1000倍の計算量が必要（1000倍時間がかかる！）となります．

つまり，陽解法では，空間解像度を細かくしたい場合に，時間刻み幅もそれにあわせて変更する必要があります．
できることならそういう事を気にせずに（全く気にしなくても良いわけではありませんが・・）自由に空間解像度を上げ下げしたいものです．
その要求を満たす解法として陰解法と呼ばれる方法があり，陰解法は無条件安定であることがわかっています（安定な計算に対するαに関する条件が無い）．
（演習では時間の関係から陰解法についてはできませんが，各人陰解法のプログラムの作成をしてみることをおすすめします．陰解法では連立一次方程式を解く必要がありますが，それらのプログラムは良いプログラミングの練習にもなります．）

レポート課題

次のノイマン条件下での拡散反応方程式について、

(2) $\begin{equation*} \left\{ \begin{array}{lc} \frac{\partial u(x, t)}{\partial t} = D\frac{\partial^2u(x, t)}{\partial x^2} & (x \in (0, \pi), 0.5>t> 0) \\ u(x, 0) = \cos(x) + \cos(7x) & ( x \in [0, \pi] ) \\ \frac{\partial u(0, t)}{\partial x} = A, \frac{\partial u(\pi, t)}{\partial x} = A \end{array} \right. \end{equation*}$

拡散係数 $D$ を $D=0.1, 1.0, 10.0$ と三つの異なる場合について、同時にアニメーションに描画する（t=0から0.5までのアニメーション，A=0.001として，時刻tもwindow内に表示）プログラムrep20190522.py(もしくはipynb)作成せよ。
異なる拡散係数を同時に表示するので、安定性条件に留意してすべての場合で計算が破綻しないように適切なパラメータを選択してください。
また、色の違いで拡散係数の違う拡散方程式の解を表示し、その情報も載せること。

採点基準はこれまで同様、コードの可読性、アニメーションの表示方法、必要要件の有無が採点対象となります。

提出期限：2019年5月30日12:00まで

おまけ１（反応拡散系の初歩）

(1) 初めに出てきた偏微分方程式をもう一度考えます.

(3) $\begin{equation*} \left\{ \begin{array}{lc} \frac{\partial u(x, t)}{\partial t} = D\frac{\partial^2u(x, t)}{\partial t^2} + f(x, t) & (0 < x < \ell, t> 0) \\ u(0, t) = a, u(\ell, t) = b & (t>0) \\ u(x, 0) = u_0(x) & (0<x<\ell) \\ u_0(x) = a, u_0(\ell) = b & \end{array} \right. \end{equation*}$

この式を

$f(x, t) = 2x - x^3$

$L=1, D=0.1$

境界条件: $a = b = 0$

初期条件: $x \le P$ ならば $u_0(x) = 1.0, x > P$ ならば $u_0(x) = -1.0$

等とし、その様子を見よ（P をいろいろ変えてみよ）

(この式はGinzburg-Landau方程式と呼ばれるものです.)