10.2 varf, matrix を用いて、 連立1次方程式を作って解く

solve, problem を使うと、 連立1次方程式 $ Au=f$ が「見えない」けれど、以下のようにして varf を用いると係数行列 $ A$ 、右辺の既知ベクトル $ f$ が得られる。 そうするのがお勧めなのだそうだ。

(おまけ: この辺の説明を読んでいて、 ようやく固有値問題を解くコードが理解出来るようになった。)

poisson3.edp

// poisson3.edp
//   Freefem++ poisson3.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;

varf a(u,v)= int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) + on(C,u=0);
matrix A=a(Vh,Vh);

varf L(unused,v) = int2d(Th)(f*v) + on(C,unused=0);
Vh F; F[] = L(0,Vh);

u[]=A^-1*F[];

plot(u);

このプログラムでは FVh の元としたが、 次のように配列にしても良い。
real [int] F=L(0,Vh);
u[] = A^-1 * F;

連立1次方程式のソルバーとしては、次のようなものがある、 とのことであるが、普通は UMFPACK を使うのだそうだ。 どういうように使い分けるのかは、良く知らない (個々の用語は良く使われるものなので、一応は意味が分かるのだけれど、 実際にどういう条件の下でどちらを選ぶかの判断が出来ない)。
solver=   対象とする行列
UMFPACK multi-frontal LU
CG CG法 疎, 正値対称
LU   skyline, 非対称
Crout Crout法 skyline, 対称
Cholesky Cholesky法 skyline, 対称, 正値
GMRES GMRES法
sparsesolver  


反復法では、eps=$ \eps$ で停止則を指定する。 $ \eps>0$ のときは

$\displaystyle \left\Vert Ax-b\right\Vert<\frac{\eps}{\left\Vert A x_0-b\right\Vert}
$

$ \eps<0$ のときは

$\displaystyle \left\Vert Ax-b\right\Vert<\left\vert\eps\right\vert.
$

桂田 祐史
2018-01-17