C-N_matlab_CrankNicolson_

%实验八 %六点对称格式求解抛物型方程 function Err=C_N(h,tau) T=10; l=1; a=1; r=a*tau/(h*h); x=0:h:l; t=0:tau:T; J=l/h; N=T/tau; u=zeros(N+1,J+1); %% 边值定义为0 u(1,:)=sin(pi*x); %% 初值条件 %%C_N格式 row=[1:J-1,1:J-2,2:J-1]; col=[1:J-1,2:J-1,1:J-2]; val=[(1+r)*ones(1,J-1),-r/2*ones(1,J-2),-r/2*ones(1,J-2)]; A=sparse(row,col,val); row2=[1:J-1,1:J-2,2:J-1]; col2=[1:J-1,2:J-1,1:J-2]; val2=[(1-r)*ones(1,J-1),r/2*ones(1,J-2),r/2*ones(1,J-2)]; B=sparse(row2,col2,val2); for n=1:N u(n+1,2:end-1)=A\(B*u(n,2:end-1)'); end [X,T]=meshgrid(x,t); %subplot(1,2,1) %mesh(X,T,u)%数值解图像 %subplot(1,2,2) %mesh(X,T,uu(X,T))%真解图像 Eh=uu(X,T)-u;%误差 %mesh(X,T,Eh)%误差函数图像 u1=uu(0.5,1/2);%t=0.5,x=0.1时刻的真解 u2=u(end,0.5*J+1);%t=0.5,x=0.1时刻的数值解 eh=abs(u1-u2);%t=0.5,x=0.1时刻数值解与真解的误差 err=max(abs(Eh(end,:)));%t=1/2时的最大误差||e_h^N||_C errs=max((abs(Eh'))); %% 不同时间层上的误差范数 plot(t,errs) Err=max(max((abs(Eh'))));%||e_τ||_C function y=uu(x,t) y=exp(-pi*pi*t).*sin(pi*x);