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 では それほど破滅的なことにはならない。

桂田 祐史
2017-06-19