tic; clear clc M=[10,20 40 80];%空间步数 N=2*M;%时间步数 for k=1:length(M) h=1/M(k);%空间步长 tau=1/N(k);%时间步长 s=tau/h;%步长比 x=0:h:1; t=0:tau:1; y=inline('exp(x+t)','x','t');%真解函数 for i=1:length(x) for j=1:length(t) exact(i,j)=y(x(i),t(j));%真解 end end u=zeros(M(k)+1,N(k)+1);%数值解内存单元 for i=2:M(k) u(i,1)=exp(x(i));%初值u(i,0) u(i,2)=exp(x(i))+tau*exp(x(i))+(tau^2/2)*exp(x(i));%第二层值u(i,1) end u(1,:)=exp(t);%边值u(0,t) u(M(k)+1,:)=exp(1+t);%边值u(1,t) phi=zeros(M(k)-1,N(k)); for i=1:N(k) phi(1,i)=exp(t(i)); phi(M(k)-1,i)=exp(1+t(i)); end A=diag((2*(1-s^2))*ones(M(k)-1,1))+diag(s^2*ones(M(k)-2,1),1)+diag(s^2*ones(M(k)-2,1),-1); for j=2:N(k) u(2:M(k),j+1)=A*u(2:M(k),j)-u(2:M(k),j-1)+s^2*phi(:,j);%数值解 end error=abs(u(2:M(k),2:end)-exact(2:M(k),2:end));%误差 error_inf(k)=max(max(error)); %error=exact-u; subplot(1,2,1) [X,Y]=meshgrid(t(end:-1:2),x(2: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=1:length(M)-1 H=error_inf(p)/error_inf(p+1); NORM(p)=log2(H); end subplot(1,2,2) plot(1:length(M)-1,NORM,'-bp') ylabel('误差阶数'); text(2,2,'这就是误差阶数!!!!!') grid on axis square toc;
效果图: