B..3 コピペして試す

$ (i,j)$成分を指定して行列を作る、という例を一つくらい。

Hilbert行列 myhilbert.m
% うっかり作ってしまったけれど、hilb(n) という関数が標準で用意されている。
function a=myhilbert(n)
  a=zeros(n,n);
  for i=1:n
      for j=1:n
          a(i,j)=1/(i+j-1);
      end
  end
end
5次のHilbert行列の逆行列と条件数を求めてみる。
edit myhilbert
(上のコードを入力して保存した後)
h=myhilbert(10)
inv(h)
cond(h)
(逆行列の成分はほぼ $ 10^{12}$ 程度の大きさになり、 条件数もほぼ $ 1.6\times 10^{13}$ 程度になる。 -- 非常に大きい、つまり悪条件。)


1次元ラプラシアン、逆行列や LU 分解を求めてみる。

1次元ラプラシアン (Dirichlet境界条件)
L=1;
N=5;
h=L/N;
e=ones(N-1,1);
a=spdiags([-e 2*e -e],-1:1,N-1,N-1)/(h*h);
full(a)

b=inv(a);
full(b)

[l u]=lu(a);
full(l)
full(u)

2次元ラプラシアン (Dirichlet境界条件)
W=2;
H=1;
Nx=8;
Ny=4;
hx=W/Nx;
hy=H/Ny;
m=Nx-1;
n=Ny-1;
ex=ones(m,1);
ax=spdiags([-ex 2*ex -ex],-1:1,m,m)/(hx*hx);
ey=ones(n,1);
ay=spdiags([-ey 2*ey -ey],-1:1,n,n)/(hy*hy);
Im=speye(m,m);
In=speye(n,n);
a=kron(ay,Im)+kron(In,ax);
full(a)

ans =

    64   -16     0     0     0     0     0   -16     0     0     0     0     0     0     0     0     0     0     0     0     0
   -16    64   -16     0     0     0     0     0   -16     0     0     0     0     0     0     0     0     0     0     0     0
     0   -16    64   -16     0     0     0     0     0   -16     0     0     0     0     0     0     0     0     0     0     0
     0     0   -16    64   -16     0     0     0     0     0   -16     0     0     0     0     0     0     0     0     0     0
     0     0     0   -16    64   -16     0     0     0     0     0   -16     0     0     0     0     0     0     0     0     0
     0     0     0     0   -16    64   -16     0     0     0     0     0   -16     0     0     0     0     0     0     0     0
     0     0     0     0     0   -16    64     0     0     0     0     0     0   -16     0     0     0     0     0     0     0
   -16     0     0     0     0     0     0    64   -16     0     0     0     0     0   -16     0     0     0     0     0     0
     0   -16     0     0     0     0     0   -16    64   -16     0     0     0     0     0   -16     0     0     0     0     0
     0     0   -16     0     0     0     0     0   -16    64   -16     0     0     0     0     0   -16     0     0     0     0
     0     0     0   -16     0     0     0     0     0   -16    64   -16     0     0     0     0     0   -16     0     0     0
     0     0     0     0   -16     0     0     0     0     0   -16    64   -16     0     0     0     0     0   -16     0     0
     0     0     0     0     0   -16     0     0     0     0     0   -16    64   -16     0     0     0     0     0   -16     0
     0     0     0     0     0     0   -16     0     0     0     0     0   -16    64     0     0     0     0     0     0   -16
     0     0     0     0     0     0     0   -16     0     0     0     0     0     0    64   -16     0     0     0     0     0
     0     0     0     0     0     0     0     0   -16     0     0     0     0     0   -16    64   -16     0     0     0     0
     0     0     0     0     0     0     0     0     0   -16     0     0     0     0     0   -16    64   -16     0     0     0
     0     0     0     0     0     0     0     0     0     0   -16     0     0     0     0     0   -16    64   -16     0     0
     0     0     0     0     0     0     0     0     0     0     0   -16     0     0     0     0     0   -16    64   -16     0
     0     0     0     0     0     0     0     0     0     0     0     0   -16     0     0     0     0     0   -16    64   -16
     0     0     0     0     0     0     0     0     0     0     0     0     0   -16     0     0     0     0     0   -16    64

>>



桂田 祐史