まずは 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 |
-->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 |
計算時間 (秒) | |
10 | 0.05 |
20 | 2.17 |
25 | 8.45 |
30 | 28.69 |
40 | 191.82 |