tic;

clear

clc

M=[,  ];%空间步数

N=*M;%时间步数

for k=:length(M)

h=/M(k);%空间步长

tau=/N(k);%时间步长

s=tau/h;%步长比

x=:h:;

t=:tau:;

y=inline('exp(x+t)','x','t');%真解函数

for i=:length(x)

    for j=:length(t)

        exact(i,j)=y(x(i),t(j));%真解

    end

end

u=zeros(M(k)+,N(k)+);%数值解内存单元

for i=:M(k)

    u(i,)=exp(x(i));%初值u(i,)

    u(i,)=exp(x(i))+tau*exp(x(i))+(tau^/)*exp(x(i));%第二层值u(i,)

end

u(,:)=exp(t);%边值u(,t)

u(M(k)+,:)=exp(+t);%边值u(,t)

phi=zeros(M(k)-,N(k));

for i=:N(k)

    phi(,i)=exp(t(i));

    phi(M(k)-,i)=exp(+t(i));

end

A=diag((*(-s^))*ones(M(k)-,))+diag(s^*ones(M(k)-,),)+diag(s^*ones(M(k)-,),-);

for j=:N(k)

    u(:M(k),j+)=A*u(:M(k),j)-u(:M(k),j-)+s^*phi(:,j);%数值解

end

error=abs(u(:M(k),:end)-exact(:M(k),:end));%误差

error_inf(k)=max(max(error));

%error=exact-u;

subplot(,,)

[X,Y]=meshgrid(t(end:-:),x(:M(k)));

mesh(X,Y,error);

xlabel('t');

ylabel('x');

zlabel('error');

grid on

% pan on;

% zoom on;

pause(0.05)

hold on

end

legend('h=1/10,tau=1/20','h=1/20,tau=1/40','h=1/40,tau=1/80','h=1/80,tau=1/160');

title('Numerical Result')

for p=:length(M)-

H=error_inf(p)/error_inf(p+);

NORM(p)=log2(H);

end

subplot(,,)

plot(:length(M)-,NORM,'-bp')

ylabel('误差阶数');

text(,,'这就是误差阶数!!!!!')

grid on

axis square

toc;