2.2 行列の等比級数

今度は正方行列 $A$ (ただし $A$ のスペクトル半径 は $1$ より小さいとする) の等比級数

$$I+A+A^2+A^3+\cdots =\sum_{n=0}^\infty A^n$$   $$\mbox{($I$ は単位行列)}$$

を考える。 この級数は解析学では ノイマンNeumann 級数と呼ばれる。

以下の実験は上の続きで行うことを想定している。 MATLAB を終了してしまっている場合は、

この節の実験に先立ち必要なこと

a=rand(5,5)
r=max(abs(eig(a)))
a=a/r
a=0.9*a

を実行しておこう。

Neumann 級数の部分和

$$S_n=\sum_{k=0}^n A^k$$

を計算する MATLAB プログラム neuamnn.m を用意した。

neumann.m

% neumann.m --- 行列の Neumann 級数 (等比級数) の第 N 部分和
function s = neumann(a,N)
  [m,n] = size(a);
  if m ~= n
    disp('aが正方行列でない！');
    return
  end
  % 第 0 項  S_0 = I
  s = eye(n,n);
  % 第 1 項  S_1 = I + a
  t = a; s = s + t;
  % 第 2〜N 項まで加える (t が a^n になるようにしてある)
  for k=2:N
    t = t * a;
    s = s + t;
  end

» edit neumann

とすると編集ウィンドウが現れるので、 上の neumann.m の内容をコピー＆ペーストして保存して下さい (ドキュメント ライブラリの MATLAB フォルダに保存されるはずです)。 その後 neumann() が使えるようになるはずです。

» c=eye(5,5)-a
» b=neumann(a,100)
» b*c
» b=neumann(a,1000)	← もう少し精度を上げてみる
» b*c
» format long	← お望みなら表示桁数を上げて
» b*c

種明かしは次のページに。