4.1.1 陽解法

% heat1d_e.m -- 空間1次元熱方程式, 同次Dirichlet境界条件
% curl -O https://m-katsurada.sakura.ne.jp/program/fdm/heat2d.m
% https://m-katsurada.sakura.ne.jp/labo/text/new-intro-matlab.pdf

% 区間
a=0
b=1
% N等分
N=100
% N等分点
x=linspace(a,b,N+1);
% 初期値 (実は u=min(x,1-x); ですむ)
for i=0:N
    if x(i+1)<0.5
        u(i+1)=x(i+1);
    else
        u(i+1)=1-x(i+1);
    end
end
% 初期値のグラフを描いて1秒待つ
%fig=figure % 新しく図のウィンドウを出すかどうか
plot(x,u)
fig=gcf; figure(fig) % こうすると visible になる
axis([0 1 -0.1 1.1])
title('heat equation')
pause(1)

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