12.1 solve, problem

一度解けば良い定常問題などは solve() で 弱形式を定義するのと同時に解いてしまえば良い。

poisson1.edp -- solve で定義すると同時に解いてしまう

// poisson1.edp
//   Freefem++ poisson1.edp
border C(t=0,2*pi) {x=cos(t); y=sin(t);}
mesh Th = buildmesh(C(50));
fespace Vh(Th,P1);
Vh u,v;
func f=x*y;
solve Poisson(u,v,solver=LU) =	
  int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) - int2d(Th)(f*v) + on(C,u=0);
plot(u);

• dx(), dy() はそれぞれ $x$, $y$ での微分を表す。
  高階の微分は dxx(), dxy(), dyy() のようにする。
• int2d(Th) は、Th の領域全体の積分(重積分)を表す。 int2d()() は複数回使うことができる。 2つめのカッコ内に書く式には制限があるようで、簡単な式に分解して書く。
• int1d(Th,border) は境界のうち、 border の部分の 積分 (境界積分、今の場合は線積分) を表す。 int1d(Th,border1,border2) のように複数個の border が指定できる。 また int1d(Th,label1,label2) のように、 border の代わりに label を指定することもできる。 int1d()() も複数回使うことができる。
• +on(border,u=g) とか。 +on(label1,label2,u=gall) とか (+on(label1,u=g1)+on(label2,u=g2) と分けて書くことも可能)。 ベクトル値関数の場合は、 +on(border,u1=g1,u2=g2) のような複数の方程式を 書くこともできる。

problem で定義しておいて、後で呼び出して解くことも出来る。 文法は solve とほとんど同じである。

poisson2.edp -- problem で定義しておいて、 後から呼び出して解く

// poisson2.edp
//   Freefem++ poisson2.edp
border C(t=0,2*pi) {x=cos(t); y=sin(t);}
mesh Th = buildmesh(C(50));
fespace Vh(Th,P1);
Vh u,v;
func f=x*y;
problem Poisson(u,v,solver=LU) =	
  int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) - int2d(Th)(f*v) + on(C,u=0);
Poisson;
plot(u);

熱方程式を解く場合などは、 何度も弱形式を解く必要があるので、 problem の利用は有効である。 各ステップで解く連立1次方程式の係数行列は共通であるので、 LU分解は一度だけやっておけば良い。 次の例では , init=i によって、そのあたりを制御している (iが0のときだけ init が false になり、 そうでないとき init は true である。 init は「初期化済み」を意味するのでしょう。)。

heat.edp -- problem で定義しておいて…

// heat.edp --- 正方形領域で熱方程式 u_t=u_{xx}+u_{yy}+f を解く
int i,m=100;
real Tmax=1, tau=0.01, t;
func f=1;
func g1=0;
func g2=0;
func u0=sin(pi*x)*sin(pi*y);
mesh Th=square(m,m);
plot(Th,wait=true);
fespace Vh(Th,P1);
// 
Vh u=u0,up,v;
problem heat(u,v,solver=UMFPACK,init=i)=
  int2d(Th)(u*v/tau+dx(u)*dx(v)+dy(u)*dy(v))
 -int2d(Th)(f*v)
 -int1d(Th,2,3)(g2*v)
 -int2d(Th)(up*v/tau)
 +on(1,4,u=g1);
for (i=0; i < Tmax/tau; i++) {
  up=u;
  t=(i+1)*tau;
  heat;
  // plot(u,cmm="t="+t,wait=0,ps="heat"+i+".ps");
  plot(u,cmm="t="+t,wait=0);
}
plot(u,cmm="t="+t,wait=1,ps="heat.ps");