成分を指定して行列を作る、という例を一つくらい。
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 |
edit myhilbert(上のコードを入力して保存した後) h=myhilbert(10) inv(h) cond(h) |
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 >>