10.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);

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 によって、そのあたりを制御している (i0 のときだけ 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");

桂田 祐史
2018-01-17