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

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

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

  $\displaystyle \gamma_1:(x,y)=(\cos\theta,\sin\theta)$   ( $ \theta\in[0,\pi/6]$)$\displaystyle ,$    
  $\displaystyle \gamma_2:(x,y)=(\cos\theta,\sin\theta)$   ( $ \theta\in[\pi/6,\pi]$)$\displaystyle ,$    
  $\displaystyle \gamma_3:(x,y)=(\cos\theta,\sin\theta)$   ( $ \theta\in[\pi,4\pi/3]$)$\displaystyle ,$    
  $\displaystyle \gamma_4:(x,y)=(\cos\theta,\sin\theta)$   ( $ \theta\in[4\pi/3,2\pi]$)$\displaystyle .$    

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

$\displaystyle 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方程式の境界値問題を考える。

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

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

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

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

$\displaystyle \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}$ に対して線積分

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

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

上の境界値は、

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

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

$\displaystyle \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);



Subsections
桂田 祐史