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