7.1 ある晩 Octave, Scilab で遊んでみました

まずは Octave で解いてみよう。

eigen_square.m

## 正方形領域のラプラシアンの固有値
function retval = eigen_square(n)
  h = 1/n;
  B=diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
  I=eye(n-1,n-1);
  A = - n * n * (- 4 * kron(I,I) + kron(B,I) + kron(I,B));
  retval = eig(A);
endfunction

eigen_square2.m

## 正方形領域のラプラシアンの固有値、固有関数
function [v,lambda] = eigen_square2(n)
  h = 1/n;
  B=diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
  I=eye(n-1,n-1);
  A = - n * n * (- 4 * kron(I,I) + kron(B,I) + kron(I,B));
  [v,lambda] = eig(A);
endfunction

  octave:1> eigen_square(10)

Scilab では、retval=eig(A) の代わりに retval=spec(A)、[v,lambda]=eig(A) の代わりに[v,lambda]=bdiag(A) くらいで行くか？ sparse() という命令があるそうだ… 期待に胸が膨らむ。

eigen_square3.m

# 正方形領域のラプラシアンの固有値 (Scilab version) バグあり
function retval = eigen_square3(n)
  h = 1/n;
  B=sparse(diag(ones(n-2,1),1)+diag(ones(n-2,1),-1));
  I=speye(n-1,n-1);
  A = - n * n * (- 4 * kron(I,I) + kron(B,I) + kron(I,B));
  retval = spec(A);
endfunction

おや、kron(I,I) でエラーになる。

-->A = - n * n * (- 4 * kron(I,I) + kron(B,I) + kron(I,B));
    p                           !--error   246 
impossible to overload this function for given argument type(s)
       undefined function %sp_kronm

つまり、残念ながら疎行列用の kron() はないわけだ。 また疎行列用の spec() もないらしい。 結局、Octave と本質的には変わらずに

改訂版 eigen_square3.m

# 正方形領域のラプラシアンの固有値 (Scilab version)
function retval = eigen_square3(n)
  B=diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
  I=eye(n-1,n-1);
  A = n * n * (4 * kron(I,I) - kron(B,I) - kron(I,B));
  retval = spec(A);
endfunction

とするしかない。 timer();a=eigen_square3(10);timer() として計算時間を計測した。

$n$	計算時間 (秒)
10	0.05
20	2.17
25	8.45
30	28.69
40	191.82

確かに計算時間は $n^6$ に比例している。 ちなみに事前に何もしないと $n=30$ で異常終了する。 スタックがあふれたらしい。 stacksize() で調べると、double を一単位として 1M となっていた。つまり $1000$ 次程度の正方行列が覚えられるリミットであ る。そこで stacksize(50*1000*1000) としたところ、 大きな $n$ に対しても計算できるようになった。 これなら 50M 要素 = 400MB だから、主記憶 512MB の oyabun では それほど破滅的なことにはならない。