ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

bicgstab

双共役傾斜安定化法

構文

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)
[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=bx について解きます。nn 列の係数行列 A は、正方行列でなければならず、さらに大規模でスパース行列でなければなりません。列ベクトル b の長さは n でなければなりません。A は、afun(x)A*x を返すような関数ハンドル afun になることができます。

関数のパラメーター化では、必要に応じて関数 afun および以下で説明する前提条件関数 mfun に追加のパラメーターを設定する方法が説明されています。

関数 bicgstab が収束する場合、それを伝えるメッセージが表示されます。関数 bicgstab が設定された最大反復回数に達しても収束しないか、または何らかの理由で停止した場合、警告メッセージが、相対残差 norm(b-A*x)/norm(b) と、計算が停止した反復回数と共に表示されます。

bicgstab(A,b,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 を使用し、方程式系 inv(M)*A*x = inv(M)*bx について効率的に解きます。M[] の場合、関数 bicgstab は前提条件を適用しません。M は、mfun(x)M\x を返すような関数ハンドル mfun です。

bicgstab(A,b,tol,maxit,M1,M2,x0) は、初期推定を指定します。x0[] の場合、関数 bicgstab は既定の設定のゼロ ベクトルを使用します。

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

フラグ

収束

0

関数 bicgstab は希望する許容誤差 tol に反復回数 maxit 以内で収束しました。

1

関数 bicgstab は、maxit 回の反復に対して、収束しませんでした。

2

前提条件 M は悪条件です。

3

関数 bicgstab の出力結果に変化が見られません (連続する 2 回の計算で、結果がまったく同じでした)。

4

関数 bicgstab で計算されたスカラー量の 1 つが計算を続けるには、大きすぎるかまたは小さすぎます。

flag0 でない場合は、返される解 x は、反復全体に渡り計算される最小のノルム残差内にあります。flag 出力が設定されている場合、メッセージは表示されません。

[x,flag,relres] = bicgstab(A,b,...) は、相対残差 norm(b-A*x)/norm(b) も返します。flag0 の場合は、relres <= 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) を含む残差ノルムのベクトルも返します。

行列入力と bicgstab の使用

この例はまず、A と前提条件 M1 を引数として直接与えることにより、Ax = b を解きます。

以下のコードを入力するとします。

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 2e-014.

関数ハンドルと bicgstab の使用

この例では、前述の例の行列 A に代えて、行列-ベクトル積関数 afun のハンドルを使用し、前提条件 M1 に代えてバックソルブ関数 mfun のハンドルを使用します。この例は、以下を行うファイル run_bicgstab に含まれています。

  • 第 1 引数として関数ハンドル @afun を指定し、関数 bicgstab を呼び出す。

  • run_bicgstab のすべての変数が afunmfun に適用できるように、入れ子関数として afunmfun を含みます。

以下に run_bicgstab のコードを示します。

function x1 = run_bicgstab
n = 21;
b = afun(ones(n,1));
tol = 1e-12;  
maxit = 15; 
x1 = bicgstab(@afun,b,tol,maxit,@mfun);

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

次のように入力した場合、

x1 = run_bicgstab;

MATLAB® ソフトウェアが以下のメッセージを表示します。

bicgstab converged at iteration 12.5 to a solution with relative 
residual 2e-014.

前提条件子と bicgstab の使用

この例は、前提条件の使用方法を示します。

479 行 479 列の非対称実スパース行列 west0479 を読み込みます。

load west0479;
A = west0479;

真の解がすべて 1 のベクトルになるように、b を定義します。

b = full(sum(A,2));

許容誤差と最大反復回数を設定します。

tol = 1e-12;
maxit = 20;

bicgstab を使用して、要求された許容誤差と反復回数で解を求めます。

[x0,fl0,rr0,it0,rv0] = bicgstab(A,b,tol,maxit);

fl0 は、bicgstab が要求した反復回数 20 回以内に要求した許容誤差 1e-12 に収束しなかったため、1 となります。実際、bicgstab の動作が良好でないため、初期推定 (x0 = zeros(size(A,2),1)) が最適解であり、これは it0 = 0 で示される出力です。MATLAB® は残差履歴を rv0 に格納します。

bicgstab の動作をプロットします。

semilogy(0:0.5:maxit,rv0/norm(b),'-o');
xlabel('Iteration number');
ylabel('Relative residual');

プロットでは、解が収束しないことが示されます。前提条件子を使用して結果を改善できます。

A は非対称であるため、ilu で前提条件子を作成します。

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));
Error using ilu
There is a pivot equal to zero. Consider decreasing 
the drop tolerance or consider using the 'udiag' option.

MATLAB は不完全な LU を構成できません。前提条件子では無意味な特異因子が生成されるためです。

エラー メッセージで示されるように、減少させた棄却許容誤差でもう一度試すことができます。

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = bicgstab(A,b,tol,maxit,L,U);

bicgstab が相対残差 5.9829e-014 (rr1 の値) を導き出すため、fl1 は 0 です。棄却許容誤差 1e-6 の不完全 LU 分解を前提条件とした場合、3 回目の反復 (it1 の値) で相対残差が指定の許容誤差 1e-12 より小さくなります。出力 rv1(1)norm(b)、出力 rv1(7)norm(b-A*x2) になります。これは、bicgstab が半分の反復を使用するためです。

初期推定 (反復回数が 0) から始まる各反復での相対残差をプロットして、関数 bicgstab の進行状況を確認できます。

semilogy(0:0.5:it1,rv1/norm(b),'-o');
xlabel('Iteration Number');
ylabel('Relative Residual');

参考文献

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

[2] 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.

拡張機能

R2006a より前に導入