3.3.3 FreeFem++ プログラム (その2)

速度ポテンシャル $\phi$ に対する Laplace 方程式の Neumann 境界値問題で、 $\rd\Omega$ 上で $\bm{v}$ が与えられたとき

$$\Laplacian\phi=0 \Omega$$ (31)

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

を満たす $\phi$ を求めよ、というもの。

これは、 $\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\}$ のときは、 $(x,y)\in\rd\Omega$ において $\bm{n}=\begin{pmatrix}x y \end{pmatrix}$ であるから、 $\bm{v}\equiv \begin{pmatrix} 1 2 \end{pmatrix}$ とすると

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

速度ポテンシャルを求める solveLaplace.edp

// solveLaplace.edp
border Gamma2(t=0,2*pi) { x=cos(t); y=sin(t); }
int m=40;
mesh Th = buildmesh(Gamma2(m));
plot(Th, wait=1, ps="Th.eps");
savemesh(Th,"Th.msh"); // optional
fespace Vh(Th,P1);
Vh u,v;
func f=0;
func g1=0; // no use
func g2=x+2*y;
solve Poisson(u,v) =
   int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
  -int2d(Th)(f*v)
  -int1d(Th,Gamma2)(g2*v) // +on(Gamma1,u=g1) 
    ;
plot(u,ps="equipotential.eps");

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

図 4: 一様流の等ポテンシャル線と流線

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