1.5.1 Poisson方程式

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

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

$$-\Laplacian u=f \Omega ,\quad u=0 \rd\Omega .$$ (1)

ただし

$$f(x,y)=xy$$   ( $(x,y)\in\Omega$)

とする。

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

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

plot() コマンドに wait=true があると、 そこで一時停止する。前に進めるためには、 キーを叩く。最後は 入力で終了する。

FreeFEM のプログラムには、 .edp という拡張子をつけることになっている。 フランス語で偏微分方程式を équations aux dérivées partielles というのでその頭文字を取ったのだとか。 (英語では partial differential equation を PDE と訳すことがよくある。)

version 4.15 から Markdown 形式のファイルも入力できるようになったそうである。

境界値問題 (1) の弱形式は

$$\int_\Omega\left(u_xv_x+u_yv_y\right)\DxDy=\int_{\Omega}f v\;\DxDy$$   ($v\in X$)

である ($f=0$ の Laplace 方程式の場合、 前回講義の Dirichlet 原理のところで一瞬顔を見せた式である。)。 弱形式については、今後の講義で詳しく解説する。 次のプログラムの20行目にこの式が「書いてある」。

19-21行で一つの命令であるが、何とこれで有限要素解が求まり、 以下は結果の表示をしているだけである。

テキスト・エディターで poisson.edp を編集してから、 FreeFem++ poisson.edp を再実行してみると良い。

poisson.edp

   1 // poisson.edp
   2 // https://m-katsurada.sakura.ne.jp/program/fem/poisson.edp
   3 // https://m-katsurada.sakura.ne.jp/labo/text/welcome-to-freefem/node4.html
   4 
   5 // 境界の定義 (単位円), いわゆる正の向き
   6 border Gamma(t=0,2*pi) { x=cos(t); y=sin(t); }
   7 // 三角形要素分割を生成 (境界を50に分割) 
   8 mesh Th = buildmesh(Gamma(50));
   9 plot(Th,wait=true); // plot(Th,wait=true,ps="Th.eps");
  10 // 有限要素空間は P1 (区分的１次多項式) 要素
  11 fespace Vh(Th,P1);
  12 Vh u,v;
  13 // Poisson 方程式 -△u=f の右辺
  14 func f = x*y;
  15 // 現在時刻をメモ
  16 real start = clock();
  17 // 問題を解く
  18 solve Poisson(u,v)
  19  = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v)
  20    +on(Gamma,u=0);
  21 // 可視化 (等高線)
  22 plot(u,wait=true);
  23 // 塗りつぶし
  24 real [int] levels =-0.012:0.001:0.012;
  25 //plot(u,viso=levels,fill=true,wait=true);
  26 // 可視化 (3次元) --- マウスで使って動かせる
  27 //plot(u,dim=3,wait=true);
  28 plot(u,dim=3,viso=levels,fill=true,wait=true);
  29 // 計算時間を表示
  30 cout << " CPU time= " << clock() - start << endl;

図 1: 三角形分割
図 2: 等高線と鳥瞰図