2.1 Laplace方程式を解いてみよう

次のような Laplace 方程式の境界値問題を考えよう。

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

$$\gamma_1:(x,y)=(\cos\theta,\sin\theta) \theta\in[0,\pi/6] ,$$

$$\gamma_2:(x,y)=(\cos\theta,\sin\theta) \theta\in[\pi/6,\pi] ,$$

$$\gamma_3:(x,y)=(\cos\theta,\sin\theta) \theta\in[\pi,4\pi/3] ,$$

$$\gamma_4:(x,y)=(\cos\theta,\sin\theta) \theta\in[4\pi/3,2\pi] .$$

$v_n\colon\rd\Omega\to\mathbb{R}$ を

$$v_n:= \left\{ \begin{array}{ll} -2 & \text{(on $\gamma_1$)} ... ...n $\gamma_3$)}\\ 0 & \text{(on $\gamma_2\cup\gamma_4$)} \end{array} \right.$$

で与える。

次のLaplace方程式の境界値問題を考える。

$$\Laplacian u=0$$   in $\Omega$$$,\quad \frac{\rd u}{\rd n}=v_n$$   on $\rd\Omega$

これは実は2次元非圧縮ポテンシャル流の速度ポテンシャル $u$ を求める問題である。 速度場を $\bm{v}$ とすると、境界上で

$$\frac{\rd u}{\rd n}=\bm{v}\cdot\bm{n}$$

が成り立つことが知られているので、

$$\bm{v}\cdot\bm{n}= \left\{ \begin{array}{ll} -2 & \text{(on $\... ...n $\gamma_3$)}\\ 0 & \text{(on $\gamma_2\cup\gamma_4$)} \end{array} \right.$$

が成り立っている、ということである。 曲線$C$と法線ベクトル$\bm{n}$ に対して線積分

$$\int_{C}\bm{v}\cdot\bm{n}\;\D s$$

は、単位時間に $C$ を超えて流れる流体の量 (面積) を意味する。

上の境界値は、

$$\int_{\rd\Omega}\frac{\rd u}{\rd n}\;\D s=0 \heartsuit$$

を満たすように選んでいる。実はこの条件は、 境界値問題の解が存在するための必要十分条件である。 実際、必要性は次のようにして分かる:

$$\int_{\rd\Omega}\frac{\rd u}{\rd n}\;\D s =\int_{\rd\Omega}\bm{v... ...\bm{n}\;\D s =\int_{\Omega}\Div\bm{v}\;\D\bm{x} =\int_{\Omega}0\;\D\bm{x}=0.$$

この問題を解くプログラムを1から作成しても良いが、 Laplace方程式やPoisson方程式のプログラムを探して、 それを書き換えることでプログラムを作成してみよう。

(Poisson方程式 $-\Laplacian u=f$ で $f=0$ とすると、 Laplace方程式であることに注意。)

Poisson方程式用のプログラム poisson-kikuchi.edp は公開してある。それを叩き台にしてみる。

curl -O https://m-katsurada.sakura.ne.jp/program/fem/poisson-kikuchi.edp
cp poisson-kikuchi.edp my-laplace.edp

この my-laplace.edp を書き換えて行こう。

これは Poisson 方程式のDirichlet-Neumann境界値問題を解くプログラムであるが、 $f=0$ とすれば Laplace 方程式となる。

+on(Gamma1,) を削除すると Dirichlet境界条件はなくすことができる。

$\dsp\int_{\gamma_1\cup\gamma_2\cup\gamma_3\cup\gamma_4}v_n v\:\D s$ は、各パートにわけて与えればよい。例えば -int1d(Th,gamma1)(g21*v)-int1d(Th,gamma3)(g23*v) とする。$\gamma_2$ と $\gamma_4$ では $v_n=0$ なので、 積分を計算する必要はない。 $\gamma_1$, $\gamma_3$ での $g_2$ を g21, g23 で与えることにした。

解 $u$ の等高線は、いわゆる等ポテンシャル線であるが、 その gradient は速度場 $\bm{v}$ なので、 それを描くと良いかもしれない。

// ベクトル場の表示
Vh u1,u2;
u1=dx(u);
u2=dy(u);
plot([u1,u2],wait=1);

解答例

// my-laplace.edp
// https://m-katsurada.sakura.ne.jp/program/fem/poisson-kikuchi.edp を元にした
// 菊地文雄, 有限要素法概説, サイエンス社

border gamma1(t=0,pi/6)      { x=cos(t); y=sin(t);  }
border gamma2(t=pi/6,pi)     { x=cos(t); y=sin(t); }
border gamma3(t=pi,4*pi/3)   { x=cos(t); y=sin(t);  }
border gamma4(t=4*pi/3,2*pi) { x=cos(t); y=sin(t); }
int m=10;
mesh Th = buildmesh(gamma1(m)+gamma2(5*m)+gamma3(2*m)+gamma4(4*m));
plot(Th,wait=true,ps="Th.eps");
savemesh(Th,"Th.msh");
fespace Vh(Th,P1);
Vh u,v;
func f=0;
func g21=-2;
func g23=1;
solve Poisson(u,v)=
  int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
 -int2d(Th)(f*v)
 -int1d(Th,gamma1)(g21*v)
 -int1d(Th,gamma3)(g23*v);
plot(u,wait=1,ps="my-laplace.eps");

//3次元鳥瞰図
//real [int] levels =0.0:0.01:1.0;
//plot(u,dim=3,viso=levels,fill=true,wait=true);

// ベクトル場の表示
Vh u1,u2;
u1=dx(u);
u2=dy(u);
plot([u1,u2],wait=1);