Matlab:Crank Nicolson方法求解线性抛物方程

tic;
clear
clc
M=[10,20,40,80,160,320,640];%x的步数
K=M; %时间t的步数
for p=1:length(M)
hx=1/M(p);
ht=1/K(p);
r=ht/hx^2; %网格比
x=0:hx:1;
t=0:ht:1;
numerical=zeros(M(p)+1,K(p)+1);
numerical(:,1)=exp(x); %初始值
numerical(1,:)=exp(t); %边值
numerical(M(p)+1,:)=exp(t+1); %边值
a=-r/2*ones(M(p)-2,1);b=(1+r)*ones(M(p)-1,1);c=-r/2*ones(M(p)-2,1);
fun1=inline('exp(x+t)','x','t');
for i=1:length(x)
    for j=1:length(t)
Accurate(i,j)=fun1(x(i),t(j));
    end
end
d=r/2*ones(M(p)-2,1);e=(1-r)*ones(M(p)-1,1);f=r/2*ones(M(p)-2,1);
B=diag(d,-1)+diag(e,0)+diag(f,1);
fun2=inline('0','x','t');
for i=1:M(p)-1
for k=1:K(p)
    f(i,k)=ht*fun2(x(i+1),t(k)+ht/2);
end
end
for k=1:K(p)
f(1,k)=r/2*(numerical(1,k+1)+numerical(1,k));
f(M(p)-1,k)=r/2*(numerical(M(p)+1,k+1)+numerical(M(p)+1,k));
end
for k=1:K(p)
    right_vector=f(:,k)+B*numerical(2:M(p),k);
    numerical(2:M(p),k+1)=chase(a,b,c,right_vector);
end
error=numerical(2:M(p),2:K(p))'-Accurate(2:M(p),2:K(p))';
error_inf(p)=max(max(error));
figure(p)
[X,Y]=meshgrid(x,t);
subplot(1,3,1)
mesh(X,Y,Accurate');
xlabel('x'),ylabel('t');zlabel('Accurate');
title('the image of Accurate result');
grid on
subplot(1,3,2)
mesh(X,Y,numerical');
xlabel('x'),ylabel('t');zlabel('numerical');
title('the image of numerical result');
grid on
subplot(1,3,3)
mesh(X,Y,numerical'-Accurate');
xlabel('x'),ylabel('t');zlabel('error');
title('the image of error result');
grid on
end
for k=2:length(M)
    H=error_inf(p-1)/error_inf(p);
E_inf(k-1)=log2(H);
end
figure(length(M)+1)
plot(1:length(M)-1,E_inf,'-r v');
ylabel('误差阶数');
title('Crank nicolson 误差阶数');
grid on
toc;