9.2 陰解法

% heat1d_i.m -- 空間1次元熱方程式, 同次Dirichlet境界条件

% 区間
a=0;
b=1;
% N等分
N=100;
% N等分点
x=linspace(a,b,N+1);
% 初期値
u=min(x,1-x);
% 初期値のグラフを描いて1秒待つ
%fig=figure
plot(x,u);
fig=gcf; figure(fig) % こうすると visible になる
axis([0 1 -0.1 1.1]);
title('heat equation');
pause(1);

tMax=1;
h=(b-a)/N;
lambda=0.5;
tau=lambda*h*h;
theta=0.5; % Crank-Nicolson
a=sparse((1+2*theta*lambda)*eye(N-1,N-1)-theta*lambda*(diag(ones(N-2,1),1)+diag(ones(N-2,1),-1)));

nMax=tMax/tau;
% t が dt 増えるごとに描画
dt = 0.01;
nskip = round(dt / tau);
for n=1:nMax
    f=(1-2*(1-theta)*lambda)*u(2:N)+(1-theta)*lambda*(u(1:N-1)+u(3:N+1));
    u(2:N)=a\f';
    u(1)=0;
    u(N+1)=0;
    if mod(n,nskip)==0
      plot(x,u);
      axis([0 1 -0.1 1.1]);
      title(['heat equation, t=', num2str(tau*n, '%4.2f')]) 
      pause(0.1);
    end
end