% % Crank-Nicolson method for the heat equation % with homogeneous Dirichlet boundary conditions % elementary code which do not exploit Matlab peculiarities % % Shows the order of the error for this method in time (wrt k, the time step) % clear all; global A mm = 10; h = 1/mm ; m=mm-1; x = [h:h:1-h] ; n = 10 ; %interval of time [0,T]=[0,1]; T=1; k = T/n ; u = sin(pi*x)' ; % u = ones(size(x))' ; U = [0;u;0]; r=k/(2*h^2); e=ones(m,1); A=1/(h^2)*spdiags([e -2*e e],-1:1,m,m); options=odeset('AbsTol',1e-10,'RelTol',1e-10); [tt uu]=ode15s('Au',[0 1],u,options); nexp=5; kk=[]; for j=1:nexp n=n*2; u = sin(pi*x)' ; U = [0;u;0]; k=1/n; kk=[kk k]; Aleft=speye(m)-k/2*A; % time stepping for l=1:n, rhs = u+k/2*A*u ; u = Aleft\rhs ; U = [U,[0;u;0]] ; end err(j)=norm(uu(end,:)'-U(2:end-1,end)); end figure(3) loglog(kk, err,'-o') hold on loglog(kk,kk.^2,'k') title('loglog plot off the norm 2 of the error vs the time step') xlabel('k') ylabel('norm of the error') hold off figure(1) t = [0:k:1] ; mesh([0 x 1],t,U') view(140, 30) set(gca,'FontSize',15) xlabel('space x') ylabel('time t') zlabel('u') figure(2) contourf([0 x 1],t,U',20) set(gca,'FontSize',15) colorbar; xlabel('space x') ylabel('time t')