3.3 実験プログラムたたき台

実験用行列生成のヒント

% ランダム Hessenberg 行列
function ret = rand_h(n)
  ret = triu(rand(n,n)) + diag(rand(n-1,1),-1);

% ランダム三重対角行列
function ret = rand_t(n)
  ret=diag(rand(n,1),0)+diag(rand(n-1,1),1)+diag(rand(n-1,1),-1);

% ランダム実対称三重対角行列
function ret = rand_st(n)
  u=rand(n-1,1);
  ret=diag(rand(n,1),0)+diag(u,1)+diag(u,-1);

% ランダム実対称行列
function ret = rand_s(n)
  a=rand(n,n);
  ret = (a+a')/2;

% QR 変換 (回数指定)
function QR = qr_iteration(a,n)
  QR = a;
  for i=1:n
    [q r]=qr(QR);
    QR=r*q;
  end

% 下三角部分のノルム
function n = norm_l(a)
  n = norm(a-triu(a));

beki.m

   1 %
   2 % beki.m --- 冪乗法の系統の算法の実験
   3 %
   4 
   5 format long
   6 
   7 a=[1,1,1;1,2,2;1,2,3]
   8 eigen=eig(a)
   9 
  10 maxiter=20; r=zeros(maxiter,1);
  11 
  12 disp("冪乗法")
  13 x=ones(3,1);
  14 for k=1:maxiter
  15   y=a*x; x=y/norm(y); r(k)=x'*(a*x);
  16 end
  17 subplot(1,4,1)
  18 semilogy(1:maxiter,abs(r-eigen(3)))
  19 % 絶対値最大の固有値に属する固有ベクトルを保存しておく
  20 u=x;
  21 
  22 disp("逆反復法")
  23 x=ones(3,1);
  24 for k=1:maxiter
  25   y=a\x; x=y/norm(y); r(k)=x'*(a*x);
  26 end
  27 subplot(1,4,2)
  28 semilogy(1:maxiter,abs(r-eigen(1)))
  29 
  30 disp("シフト法")
  31 kinji=0.6;
  32 ap=a-kinji*eye(3,3);
  33 x=ones(3,1);
  34 for k=1:maxiter
  35   y=ap\x; x=y/norm(y); r(k)=x'*(ap*x)+kinji;
  36 end
  37 subplot(1,4,3)
  38 semilogy(1:maxiter,abs(r-eigen(2)))
  39 
  40 x=ones(3,1);
  41 x=x-(u'*x)*u;
  42 for k=1:maxiter
  43   y=a*x; x=y/norm(y); x=x-(u'*x)*u; r(k)=x'*(a*x);
  44 end
  45 subplot(1,4,4)
  46 semilogy(1:maxiter,abs(r-eigen(2)))

householder.m

   1 % Householder
   2 function a=householder(A)
   3   [n,n] = size(A);
   4   a = A;
   5   for k=1:n-1
   6     b = a(k+1:n,k);
   7     B = a(k+1:n,k+1:n);
   8     s = - sign(b(1)) * norm(b);
   9     w = b;
  10     w(1) = w(1) - s;
  11     normw2 = 2 * s * (s - b(1));
  12     p = (B * w) / (normw2/2);
  13     alpha = (w'*p) / normw2;
  14     q = p - alpha * w;
  15     a(k+1:n,k+1:n) = B - w * q' - q * w';
  16     a(k+2:n,k) = zeros(n-k-1,1);
  17     a(k,k+2:n) = zeros(n-k-1,1)';
  18     a(k+1,k) = s;
  19     a(k,k+1) = s;
  20   end

なお、授業のレポートである福澤 [7] も参考になるかも。