4.4 (1), (2) を解くための FreeFem++ のプログラム

速度ポテンシャルを求める境界値問題

$$\Laplacian\phi=0 \Omega$$ (再1)

$$\frac{\rd\phi}{\rd\bm{n}}=\bm{v}\cdot\bm{n} \rd\Omega$$ (再2)

の場合は、 $\Gamma_1=\emptyset$ , $\Gamma_2=\rd\Omega$ , $f=0$ , $g_2=\bm{v}\cdot\bm{n}$ .

特に $\Omega=\left\{(x,y)\in\mathbb{R}^2\mid x^2+y^2<1\right\}$ , $\bm{v}\equiv \begin{pmatrix} 1 2 \end{pmatrix}$ のときは、 $(x,y)\in\rd\Omega$ において $\bm{n}=\begin{pmatrix}x y \end{pmatrix}$ であるから、

$$g_2=\bm{v}\cdot\bm{n}= \begin{pmatrix} 1 2 \end{pmatrix} \cdot \begin{pmatrix} x y \end{pmatrix} =x+2y.$$

速度ポテンシャルを求める potential2d-v0.edp ← これを保存する

// potential2d-v0.edp
//   http://nalab.mind.meiji.ac.jp/~mk/complex2/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");

fespace Vh(Th,P1);
Vh u, v, phi, psi;
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);

// ベクトル場 (u,v)=∇φ を描く (ちょっと雑なやり方)
u=dx(phi);
v=dy(phi);
plot([u,v],ps="vectorfield.eps",wait=1);

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

プログラムはテキスト・エディター (mi, emacs, Xcode, テキスト・エディット²など) で作成し、ターミナルから、

こんなふうにして実行

FreeFem++ potential2d.edp

とタイプして実行できる。リターン・キーを打つごとに次の図に進み、 最後 (ベクトル場と等ポテンシャル線を描いてお終い) は Escape キーを打って終了する。

図 6: $\Omega$ の三角形分割

図 7: 一様流の等ポテンシャル線, 流線, ベクトル場

このプログラムは、 図の PostScript データ ``counterpotential.eps'', ``vectorfield.eps'', ``both.eps'' を出力している。

(実は、上の弱形式は解の一意性がないので、 このプログラムは少し危ういところがある。)

FreeFem++ はポテンシャル問題に限らず、色々な問題を解くことが出来る。 時間が余った時の脇道用: Navier-Stokes 方程式による、 有名な円柱の周りの水の流れのシミュレーション・プログラム NS-cylinder.edp (鈴木 厚先生の Finite element programming by FreeFem++ - intermediate course から)