5 サンプルプログラム potential2d-2026.edp 解説

potential2d-2026.edp

   1 // potential2d-2026-naive.edp --- 2次元非圧縮ポテンシャル流
   2 //   https://nalab.mind.meiji.ac.jp/~mk/complex2/potential2d-2026-naive.edp
   3 //     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く
   4 
   5 
   6 border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
   7 int m=40;
   8 mesh Th=buildmesh(Gamma(m));
   9 plot(Th, wait=1, ps="Th.eps");
  10 // 次の2行は区分1次多項式を使うという意味
  11 fespace Vh(Th,P1);
  12 Vh phi, v;
  13 // 境界条件の設定
  14 //func Vn=x+2*y; // これは一様流 v=(1,2) の場合、v・n=(1,2)・(x,y)=x+2y
  15 func Vn=((x>0&&y>0)||(x<0&&y<0))*(x+2*y);//右上と左下のみ出入りがある
  16 // 整合条件 ∫ V_n ds=0 のチェック
  17 cout << "check the compatibility condition: " << int1d(Th,Gamma)(Vn) << endl;
  18 
  19 // 有限要素法で速度ポテンシャルφを求める
  20 solve Laplace(phi,v,solver=CG) =
  21   int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn*v);
  22 
  23 // 平均を0にする (解の一意性がないので、時々生じるずれを消す)
  24 real meanphi = int2d(Th)(phi)/Th.area; // Th.areaはint2d(Th)(1);でも計算可能
  25 cout << "mean phi=" << meanphi << endl;
  26 phi=phi-meanphi;
  27 
  28 // 速度ポテンシャルφの等高線 (等ポテンシャル線) を描く
  29 plot(phi,ps="contourpotential.eps",wait=1);
  30 
  31 // ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
  32 Vh v1, v2;
  33 v1=dx(phi); v2=dy(phi);
  34 plot([v1,v2],ps="vectorfield.eps",wait=1);
  35 
  36 // 等ポテンシャル線とベクトル場を同時に描く
  37 plot([v1,v2],phi,ps="both.eps", wait=1);

図 1: potential-2d-2026.edp の結果 (左下と右上で出入りがある)

• 6行目で領域$\Omega$の境界 $\Gamma:=\rd\Omega$を指定している (それで $\Omega$ が決まる)。
• 8行目で、$\Gamma$を$m$分割して、$\Gamma$の囲む範囲を三角形分割して、 それを mesh 型の変数 Th に代入している。
• 9行目で三角形分割の様子を画面に表示する。
• 11行目、有限要素空間 $V_h$ を区分的1次関数の空間と定義している。 この辺についてはこの講義では説明を省略する (知りたい人は有限要素法のテキストを読んで下さい)。 ここを変更する必要はない。
• 12行目、解 $\phi$ と試験関数 $v$ を Vh の要素とする。
• 15行目、境界値 $V_n$ の設定をしている。ここでは
  $$V_n(x,y)= \left\{ \begin{array}{ll} x+2y &\text{($(x>0\text{�... ...ext{かつ}y<0)$のとき)}\\ 0 & \text{(それ以外)} \end{array} \right.$$
  としている。(後の例 sample0.edp では $V_n(x,y)=x+2y =\begin{pmatrix}1 2 \end{pmatrix}\cdot\bm{n}$ としていて、 これは $\bm{v}=\begin{pmatrix}1 2 \end{pmatrix}$ という一様流が解となる。) これは確かに $\dsp\int_{\rd\Omega}V_n \;\D\sigma=0$ を満たしている。
• 17行目 $\dsp\int_{\rd\Omega}V_n\;\D\sigma$ を計算して表示する。 これが0になるか(0に十分近いか)どうか。
  このようにプログラムでチェックすることを勧める。
• 20〜21行目、弱形式を定義すると同時に有限要素解を求めている。 係数行列が特異なので、 連立1次方程式は Gauss の消去法のような直接法では解けない。 ここではCG法 (共役勾配法) を指定しているが、この辺は変更する必要はない。
• 24〜26行目、Neumann境界値問題の解はは定数だけの不定性を持つため ($\phi$ が解ならば $\phi+$定数 も解)、ごく稀に大きなズレが生じ、 解の描画などで問題が生じる。平均値 $\dfrac{\dsp\dint_\Omega\phi(x,y)\DxDy} {\dsp\dint_{\Omega}\DxDy}$ を引くことで、平均$=0$にしておく。 やらなくても動く場合が多い。
• 29行目、解$\phi$の等高線 (等ポテンシャル線) を描く。
• 32〜34行目、速度場 $\grad\phi$ を (矢印で) 表示する。
• 37行目、解$\phi$の等高線と速度場 $\grad\phi$ を一つの図に描く。