6.6 ポテンシャル問題の数値解法 (1) 有限要素法

ポテンシャル問題を数値的に解くことを考えよう。 この「応用複素関数」では、 有限要素法と 基本解の方法を簡単に紹介する。

差分法で解くこともできるが、 長方形領域でない問題を解くには工夫が必要になり、あまり便利でない。

有限要素法の主たるアイディアは次の2つ:

(1)

弱形式を用いる。

(2)

領域を三角形、四面体などの有限要素に分割し、 近似解や試験関数に区分的多項式を採用する。

この講義では有限要素法の詳細は解説できないが、 幸いFreeFEM というソフトを用いると、 弱形式さえ分かれば、 有限要素についてはソフトに任せにして、数値計算ができる。

実は Dirichlet 原理の証明中に現れた (6.8) は Laplace方程式の Dirichlet 境界値問題の弱形式である。 (弱形式については、次回解説を行う。)

今回は「百聞は一見にしかず」で、 まずはプログラム(スライド1枚)を紹介する。

2,3行書き換えるだけで「自分の問題」が解ける。

// potential2d-v0.edp ---  2次元非圧縮ポテンシャル流
//  速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く
border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
int m=40;
mesh Th=buildmesh(Gamma(m));
plot(Th, wait=1, ps="Th.eps");
// 次の2行は区分1次多項式を使うという意味
fespace Vh(Th,P1);
Vh phi, v, v1, v2;
// 境界条件の設定
func Vn=x+2*y; // Ωが単位円で, V=(1,2) のとき V・n=x+2y

// 速度ポテンシャルφを求め、その等高線 (等ポテンシャル線) を描く
solve Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn*v);
plot(phi,ps="contourpotential.eps",wait=1);
// ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
v1=dx(phi); v2=dy(phi);
plot([v1,v2],ps="vectorfield.eps",wait=1);

// 等ポテンシャル線とベクトル場を同時に描く
plot([v1,v2],phi,ps="both.eps", wait=1);