1 はじめに

(内部での熱発生を考慮した) 熱方程式の初期値境界値問題

に対して ($\Omega$, $\Gamma_1$, $\Gamma_2$ については、 Poisson方程式の例題のそれと同じとする)、 時間微分を$\theta$法で離散化したプログラムを作って、 安定性の条件を調べてみよう。

$\theta$法では、 $0\le\theta\le 1$を満たす$\theta$ を取って、 微分方程式 $\dfrac{\rd u}{\rd t}(\cdot,t_n) =\Laplacian u(\cdot,t_n)+f$ を

$$\frac{1}{\Delta t}\left(u^{n+1}-u^n\right)=\Laplacian u^{n+\theta}+f$$ (2)

で置き換えるわけである。ここで

$$u^n:=u(\cdot,t_n),\quad t_n:=n\Delta t,$$ (3)

$$u^{n+\theta}:=\theta u^{n+1}+(1-\theta)u^n=\theta u(\cdot,t_{n+1})+(1-\theta)u(\cdot,t_n).$$ (4)

(4) を (2) に代入して

$$\frac{1}{\Delta t}\left(u^{n+1}-u^n\right) =\theta\Laplacian u^{n+1}+(1-\theta)\Laplacian u^{n}+f.$$ (5)

これから弱形式を導く。

後退Euler法については、弱形式とサンプル・プログラム heatB.edp を与えておいた。 $\theta$ 法については、弱形式はあるが、プログラムがない。 このような場合にプログラムを作る作業はよくあるタイプの作業である。

heatB.edp を叩き台にして heatT.edp を作る

curl -O https://m-katsurada.sakura.ne.jp/program/fem/heatB.edp
cp heatB.edp heatT.edp