Octave のサンプル・プログラム

冪乗法

format long             % 表示桁数を多く取る
a=[1 1 1;1 2 2;1 2 3]
eigen=eig(a)            % eig() で固有値を求める (カンニング)
x=ones(3,1);            % 初期ベクトルを (1,1,1) とする
maxiter=10;
r=zeros(maxiter,1);     % 近似固有値を記録するベクトルの用意
for k=1:maxiter
  y=a*x;
  x=y/norm(y);
  r(k)=x'*(a*x);        % Rayleigh 商で近似固有値を計算し、r() に記録
end
error=r-eigen(3)        % 近似固有値の誤差を求める
plot(1:maxiter,abs(error))     % 反復回数と誤差をプロット
loglog(1:maxiter,abs(error))   % 両側対数目盛りでプロット
semilogy(1:maxiter,abs(error)) % 片側対数目盛りでプロット