1.5.2 熱伝導方程式

桂田 [4] の該当箇所
https://m-katsurada.sakura.ne.jp/labo/text/welcome-to-freefem/node11.html

領域の内部で熱が発生する (or 熱が吸収される) 場合の熱方程式の初期値境界値問題を解いてみよう。

$$\Omega:=\left\{(x,y)\in\mathbb{R}^2 \relmiddle\vert 0<x<1, 0<y<1\right\}$$

$$\Gamma_1=\left\{(0,y)\relmiddle\vert 0\le y\le 1\right\} \cup \left\{(x,0)\relmiddle\vert 0\le x\le 1\right\},$$

$$\Gamma_2=\left\{(1,y)\relmiddle\vert 0\le y\le 1\right\} \cup \left\{(x,1)\relmiddle\vert 0\le x\le 1\right\}.$$

つまり $\Omega$ は正方形領域、$\Gamma_1$ は左と下の辺を合わせたもの、 $\Gamma_2$ は右と上の辺を合わせたものである。

以下のプログラムでは

$$f(x,y)=1,\quad g_1(x,y)=0,\quad g_2(x,y)=0,\quad u_0(x,y)=\sin(\pi x)\sin(\pi y)$$

としてある。

heatB.edp

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

ターミナルでこうすれば実行できる

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

時間発展する系なので、動画で見るのが分かりやすいであろう。

https://m-katsurada.sakura.ne.jp/ana2026/heat.mp4

定常解 $v$ に収束することが分かるが、 実は、Poisson 方程式の境界値問題

$$-\Laplacian v=f$$   (in $\Omega$)$$,\quad v=g_1$$   (on $\Gamma_1$)$$,\quad \frac{\rd v}{\rd n}=g_2$$   (on $\Gamma_2$)

の解となっている (https://m-katsurada.sakura.ne.jp/labo/text/welcome-to-freefem/node8.html)。