Poisson 方程式の境界値問題で説明します (有限要素法を自習したい場合には、 菊地文雄, 『有限要素法概説』, サイエンス社, という参考書がイチオシです, 桂田研の学生で希望する人は申し出て下さい)。
内の有界な領域 (連結開集合) における Poisson 方程式の Dirichlet 問題
(PE) | ||
(DBC) |
Poisson 方程式 (PE) に、任意の ( は、 内に compact な台を持つ 級の関数全体の集合) をかけて、Green の積分公式 (多次元版部分積分) を用いると、 次式が得られます。
逆に (on ) を満たす滑らかな が、 この条件 (1) を満たせば、 は Poisson 方程式 (PE) を満たすことも容易に分かります。
微分方程式 (PE) を満たす を求める代りに、 (1) を満たす を見出すことを考えましょう。 解の存在を容易に保証できるようにするため、 関数空間を完備化するのが便利です (解を極限の論法で得よう、という解析学の常套手段を使うため)。 具体的には、 以下に定義する というソボレフ空間 (一般化された導関数を持つ関数からなる関数空間の一種) で を探すことにします。
(PE),(DBC) の弱解の定義 | |||
かつ
|
とは |
|
を三角形の「整った」合併 で近似し、 で連続で、 各三角形上で1次関数であるような関数全体を とします。 は の有限次元部分空間とみなせ、 グラフが角錐であるような関数 , , を 基底に取って、 各元をその基底の線型結合として表現することができます:
(PE),(DBC) の弱解 の定義 (書き換え) | ||||||
かつ
|
有限要素解 の定義 | ||||||
かつ
|
次のようなプログラム poisson.edp を用意します (テキスト・エディターを使って作成します。 現象数理学科学生は、CotEditor, VS code (code), Atom (atom), emacs などを使っているでしょう。 どれを使ってもOKです。)。
poisson.edp |
// poisson.edp // https://m-katsurada.sakura.ne.jp/program/fem/poisson.edp // https://m-katsurada.sakura.ne.jp/labo/text/welcome-to-freefem/node4.html // 境界の定義 (単位円), いわゆる正の向き border Gamma(t=0,2*pi) { x=cos(t); y=sin(t); } // 三角形要素分割を生成 (境界を50に分割) mesh Th = buildmesh(Gamma(50)); plot(Th,wait=true); // plot(Th,wait=true,ps="Th.eps"); // 有限要素空間は P1 (区分的1次多項式) 要素 real [int] levels =-0.012:0.001:0.012; fespace Vh(Th,P1); Vh u,v; // Poisson 方程式 -△u=f の右辺 func f = x*y; // 現在時刻をメモ real start = clock(); // 問題を解く solve Poisson(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v) +on(Gamma,u=0); // 可視化 (等高線) plot(u,wait=true); //plot(u,viso=levels,fill=true,wait=true); // 可視化 (3次元) --- マウスで使って動かせる //plot(u,dim=3,wait=true); plot(u,dim=3,viso=levels,fill=true,wait=true); // 計算時間を表示 cout << " CPU time= " << clock() - start << endl; |
ターミナルで poisson.edp を入手して実行 ($ はプロンプトなので入力しない) |
$ curl -O https://m-katsurada.sakura.ne.jp/program/fem/poisson.edp
$ FreeFem++ poisson.edp |
このプログラムを実行すると、 円盤領域の三角形分割を表示したところで停止します (plot() に wait=true としてあるため)。 先に進めるには (return) を入力します。
等高線が表示されたところでも、 先に進めるには (return) を入力します。
最後にグラフを描画しますが、 そのウィンドウを閉じるには を入力します。
キーボードからのコマンド | ||||||||||||||||||||||||||
|
plot(u,wait=true); のところを plot(u,wait=true,ps="poisson-disk.eps"); のようにすると、 画面に表示するだけでなく、PostScript 形式でファイル出力できます (TEX 文書への取り込みが簡単です)。
桂田 祐史