反復法の収束性

今回の授業では、Hilbert行列の例を利用して反復法の収束性を検討してみます。

1 準備

Hilbert行列を作成します。MATLABではhilb(n)という命令で$n\times n$のHilbert行列を作成することができます。

$$ H=( h_{ij})_{n \times n}, \quad h_{ij}=\frac{1}{j+j-1}\:. $$

n=2
A=hilb(n)

n=3
A=hilb(n)

n =  2
A =

   1.00000   0.50000
   0.50000   0.33333

n =  3
A =

   1.00000   0.50000   0.33333
   0.50000   0.33333   0.25000
   0.33333   0.25000   0.20000

Hilbert行列の性質を調べます。

• 行列式　( det(A) )
• 最大固有値と最小固有値 ( eig(A) )、さらに条件数(最大固有値/最小固有値)

n=3
A=hilb(n)
printf("行列式= %15.10f \n\n", det(A))
myeig = eig(A);
printf("固有値=\n\n")
disp(myeig)

cond_num = max(myeig)/min(myeig);

printf("\n条件数= %15.10f \n", cond_num)

n =  3
A =

   1.00000   0.50000   0.33333
   0.50000   0.33333   0.25000
   0.33333   0.25000   0.20000

行列式=    0.0004629630 

固有値=

   0.0026873
   0.1223271
   1.4083189

条件数=  524.0567775861 

2. 反復法で連立一次方程式を解く

SOR法を使って、Ax=bの解を求めてください。ここで、$b$のすべての成分を１とする。

以下のコードでは、SOR法が関数として定義されています。

特に、trilとtriuを使って、行列の下三角行列と上三角行列を取ることができます。

function [sol, errlist, iteration_num] = sor(A,b,w,max_N)

    E=tril(A,-1);
    F=triu(A,1);
    D=diag(diag(A));
    exact_sol = A\b;
        
    M=D/w+E; N=(w-1)/w*D+F;

    x=zeros(length(b),1);
    errlist = [];
    for k=1:max_N;
      x = M\(b-N*x);
      errlist(k) = norm(x - exact_sol);
      if errlist(k)<1E-5
        break
      end
    end
    sol=x;
    iteration_num =k;
end

A=hilb(n);  
b = ones(n,1);  %  Vector b can also be inputed by:  b=[1,1,1,1]';

[sol, errlist,num] = sor(A,b,1.5,1000);
sol
num
plot(log10(errlist))
grid on

sol =

    3.0000
  -24.0000
   30.0000

num =  236

課題1

• $N=3$の場合、$1.1\le w \le 1.9$の中の$w$について、近似解のL2ノルムの誤差が$10^{-6}$以下になるまでに必要な反復計算の回数を調べて、一番良い$w$を選んでください。

• $N=3,4,5$の場合、上記の計算で定めた$w$の値を利用して、SOR法で近似解を求めてください。近似解のL2ノルムの誤差が$10^{-6}$以下ことが必要です。特に、反復計算回数が20000になる可能性があるので、大きい反復計算回数で試してみてください。

• 行列の条件と反復計算の回数の関係を調べてください。

3. ガウス消去法で連立一次方程式を解く

ガウス消去法を使って、Ax=bの解を求めてください。ここで、$b$のすべての成分を１とする。

n=3
A=hilb(n)   
b = ones(n,1);  %  Vector b can also be inputed by:  b=[1,1,1,1]';

sol_x = A\b

%ここでガウス消去法によって近似解を求めてください。（ピボット選択の使用は自由です。）

n =  3
A =

   1.00000   0.50000   0.33333
   0.50000   0.33333   0.25000
   0.33333   0.25000   0.20000

sol_x =

    3.0000
  -24.0000
   30.0000

課題2

$N=3,4,5,6$の場合、Hilbert行列の条件数を調べて、近似解の誤差との関係を検討してみてください。