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

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

このプログラムでは F は Vh の元としたが、 次のように配列にしても良い。

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$ のときは

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

$\eps<0$ のときは

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