% 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 |