D. 熱方程式に対する後退 Euler 法

内部で熱が発生する (or 熱が吸収される) 場合の熱方程式の初期値境界値問題

$$u_t(x,t)=\Laplacian u(x,t)+f(x) (x,t)\in \Omega\times(0,\infty) ,$$ (11a)

$$u(x,t)=g_1(x) (x,t)\in\Gamma_1\times(0,\infty) ,$$ (11b)

$$\frac{\rd u}{\rd n}(x,t)=g_2(x) (x,t)\in\Gamma_2\times(0,\infty) ,$$ (11c)

$$u(x,0)=u_0(x) x\in\overline{\Omega}$$ (11d)

を考える。

変数 $t$ については(後退)差分近似する。すなわち $t_n:=n\Delta t$, $u^n:=u(\cdot, t_n)$ とおいて、

$$\frac{\rd u}{\rd t}(\cdot,t_n)\kinji\frac{u^n-u^{n-1}}{\Delta t}$$

という近似を採用する。

弱形式は

$$\left(\frac{u^{n}-u^{n-1}}{\Delta t},v\right) +\langle u^{n},v\rangle-(f,v)-[g_2,v]=0.$$

すなわち

$$\left(u^{n+1},v\right) -\left(u^{n},v\right) +\Delta t\langle u^{n+1},v\rangle-\Delta t(f,v)-\Delta t[g_2,v]=0.$$

既に $u^n$ が “分かっている” とき、 これは $u^{n+1}$ についての Poisson 方程式もどきに対する弱形式であり、 これまでとほぼ同様にして解くことが出来る。

heatB.edp

// heatB.edp --- 熱方程式を後退Euler法 (陰解法) で解く
// https://m-katsurada.sakura.ne.jp/program/fem/heatB.edp
// 菊地文雄, 有限要素法概説, サイエンス社のPoisson方程式の問題の非定常版
int i,m=10;
real Tmax=1, tau=0.01, t;
// 次の2行の行頭の // を削除すると、m, tau, theta が実行時に変更できる
// cout << "m dt theta: "; cin >> m >> tau >> theta;
// cout << "m=" << m << ", tau=" << tau << ", theta=" << theta << endl;
mesh Th=square(m,m);
plot(Th,wait=true);
fespace Vh(Th,P1);
func f=1; func g1=0; func g2=0;
func u0=sin(pi*x)*sin(pi*y);
Vh u=u0, uold, v;
plot(u,cmm="t=0",wait=1);
problem heat(u,v,init=i)=
  int2d(Th)(u*v+tau*(dx(u)*dx(v)+dy(u)*dy(v)))-int2d(Th)(uold*v)
 -int2d(Th)(tau*f*v)-int1d(Th,2,3)(tau*g2*v)
 +on(1,4,u=g1);
for (i=0;i<Tmax/tau;i++) {
  uold=u;
  t=(i+1)*tau;
  heat;
  plot(u,cmm="t="+t,wait=0); // ps="heat"+i+".ps" として保存することも可能
}
plot(u,wait=1);

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

ループの制御変数を i という変数にして、 problem に ,init=i と書き足すのがミソ。 最初は i が 0 であるので、init が False であるが、 それ以降は $\texttt{i}\ne 0$ であるので、init が True である、 つまり LU 分解済みだ、と指示するわけである。

時間が経過すると、 定常解に収束する(Poisson方程式の解に近く)様子が見える。