bicgstab (MATLAB Function 参考文献)

双共役傾斜安定化法

表示

x = bicgstab(A,b)
bicgstab(A,b,tol)
bicgstab(A,b,tol,maxit)
bicgstab(A,b,tol,maxit,M)
bicgstab(A,b,tol,maxit,M1,M2)
bicgstab(A,b,tol,maxit,M1,M2,x0)
bicgstab(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,...)
[x,flag] = bicgstab(A,b,...)
[x,flag,relres] = bicgstab(A,b,...)
[x,flag,relres,iter] = bicgstab(A,b,...)
[x,flag,relres,iter,resvec] = bicgstab(A,b,...)

詳細

x = bicgstab(A,b) は、線形方程式システムA*x = bをxについて解きます。n行n列の係数行列Aは、正方である必要があり、列ベクトル b は、長さがnでなければなりません。Aは、afun(x)が、A*xを出力するような関数afunです。

bicgstabが収束する場合、その影響を示すメッセージが表示されます。bicgstabが最大繰り返し回数に達した後、収束していないか、または、何らかの理由で、計算が停止した場合、ワーニングメッセージとして、相対残差 norm(b-A*x)/norm(b)を表示し、停止したときの繰り返し回数を表示します。

bicgstab(A,b,tol) は、手法の許容誤差tolを指定します。tol が、[]の場合、bicgstab はデフォルト1e-6を使います。

bicgstab(A,b,tol,maxit) は、繰り返し最大回数を指定します。maxit を []に設定すると、bicgstab は、デフォルト値 min(n,20) を使います。

bicgstab(A,b,tol,maxit,M) と bicgstab(A,b,tol,maxit,M1,M2) は、前提条件子 M、または、M = M1*M2 を使って、xに関して、システム inv(M)*A*x = inv(M)*b を効率的に解きます。 M が、[] の場合、bicgstab は、前提条件子を適用しません。M は、M\xを戻す関数です。

bicgstab(A,b,tol,maxit,M1,M2,x0) は、初期推定値を設定します。x0 が、[]の場合、bicgstab は、すべての要素がゼロのベクトルを使います。

bicgstab(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,...) は、パラメータ p1,p2,... を関数 afun(x,p1,p2,...), m1fun(x,p1,p2,...), m2fun(x,p1,p2,...) に渡します。

[x,flag] = bicgstab(A,b,...) は、収束フラグも出力します。

フラグ
収束

0
bicgstab は、繰り返し回数 maxit 以内で、希望するトレランス tol 内に収束します。

1
bicgstab は、繰り返し回数が maxit に達しましたが、収束しませんでした。

2
前提条件子 M の条件数が悪い。

3
bicgstabは、計算を進めていません(連続する2回の繰り返しの結果が同じでした)。

4
bicgstab は、関数の中に計算されたスカラ量が、繰り返しを続行するのに小さすぎる、または大きすぎるため、この手法は停止されました。

フラグ	収束
`0`	`bicgstab` は、繰り返し回数 `maxit` 以内で、希望するトレランス `tol` 内に収束します。
`1`	`bicgstab` は、繰り返し回数が `maxit` に達しましたが、収束しませんでした。
`2`	前提条件子 `M` の条件数が悪い。
`3`	`bicgstab`は、計算を進めていません(連続する2回の繰り返しの結果が同じでした)。
`4`	`bicgstab` は、関数の中に計算されたスカラ量が、繰り返しを続行するのに小さすぎる、または大きすぎるため、この手法は停止されました。

flagが0でないとき、出力された解xは、すべての繰り返しで計算された中で最小のノルム残差をもつ解です。flagの出力が指定されている場合は、メッセージは表示されません。

[x,flag,relres] = bicgstab(A,b,...) は、相対誤差norm(b-A*x)/norm(b)も出力します。flagが0の場合、これはtol以下です。

[x,flag,relres,iter] = bicgstab(A,b,...) は、xが計算された点での繰り返し回数も出力します。ここで、0 <= iter <= maxitです。iterは、整数 + 0.5で、繰り返しを通して収束の半分を示します。

[x,flag,relres,iter,resvec] = bicgstab(A,b,...) は、norm(b-A*x0)を含み、各半分の繰り返しでの残差ノルムベクトルを出力します。

例題

例題 1.

A = gallery('wilk',21);  
b = sum(A,2);
tol = 1e-12;  
maxit = 15; 
M1 = diag([10:-1:1 1 1:10]);

x = bicgstab(A,b,tol,maxit,M1,[],[]);
bicgstab converged at iteration 12.5 to a solution with relative
residual 1.2e-014

一方、つぎの行列-ベクトルの乗算の関数を使います。

function y = afun(x,n)
y = [0; 
     x(1:n-1)] + [((n-1)/2:-1:0)'; 
     (1:(n-1)/2)'] .*x + [x(2:n); 
     0];

そして、つぎの前提条件子バックソルブ関数

function y = mfun(r,n)
y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];

bicgstabへの入力として、

x1 = bicgstab(@afun,b,tol,maxit,@mfun,[],[],21);

afun と mfunは、bicgstabの外部入力n = 21を受け入れる必要があることに注意してください。

例題 2.

load west0479;
A = west0479;
b = sum(A,2);
[x,flag] = bicgstab(A,b)

bicgstabは、デフォルトの20回の繰り返しでデフォルトのトレランス 1e-6 で収束しないので、flag は 1 になります。

[L1,U1] = luinc(A,1e-5);
[x1,flag1] = bicgstab(A,b,1e-6,20,L1,U1)

上三角U1は、対角要素に一つのゼロをもつので、flag1 は 2 になります。bicgstab は最初の繰り返しで、失敗しているので、バックシュラッシュを使って、U1*y = r のようなシステムを解こうとします。

[L2,U2] = luinc(A,1e-6);
[x2,flag2,relres2,iter2,resvec2] = bicgstab(A,b,1e-15,10,L2,U2)

bicgstabは、1e-6のドロップトレランスをもち、不完全LU 分解を前提条件とした場合、6回の繰り返し(iter2の値)で、3.1757e-016のトレランス (relres2の値)に収束するので、flag2 は 0 になります。resvec2(1) = norm(b) と resvec2(13) = norm(b-A*x2)です。ユーザは、中間点と初期推定(くり返し回数 0)からスタートして、各繰り返しの最後での相対残差プロットすることで、bicgstabのプロセスをフォローすることができます。

semilogy(0:0.5:iter2,resvec2/norm(b),'-o')
xlabel('iteration number')
ylabel('relative residual')

参考

bicg, cgs, gmres, lsqr, luinc, minres, pcg, qmr, symmlq

@ (関数ハンドル), \ (バックシュラッシュ)

参考文献

Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

van der Vorst, H. A., "BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems", SIAM J. Sci. Stat. Comput., March 1992,Vol. 13, No. 2, pp. 631-644.

bicg bin2dec