微分方程数值解法编程作业二

微分方程数值解法编程作业二考虑扩散方程的初边值问题22,01,0;,0sin,01;0,1,0,0.uuxttxuxxxututt四种差分格式为：（1）向前差分格式11112nnnnjjjjuuuu（2）向后差分格式1111112nnnnjjjjuuuu（3）Crank-Nicolson格式1111111112222nnnnnnjjjjjjuuuuuu（4）DuFort-Frankel格式1111121222nnnnjjjjuuuu其中，（1）为二步显式格式，（2）、（3）为二步隐式格式，（4）为三步显式格式。（2）、（3）格式需要求解线性方程组。用四种方法计算数值解，并与初值问题的解析解相比较，得出你的结论并简要说明。解：（1）functionkuosan1%向前差分h=0.1;lambna=[0.10.5];x=1;t=0.1;J=round(x/h);x=0:h:1;u=sin(pi*x');fork=1:2N=round(t/(lambna(k)*h^2));fori=1:Nforj=2:Ju(j,i+1)=u(j,i)+lambna(k)*(u(j+1,i)-2*u(j,i)+u(j-1,i));endendu(:,N+1);w=exp(-pi^2*t)*sin(pi*x');ez=u(:,N+1)-w;table=[u(:,N+1),w,ez];disp(strcat('table',num2str(k),':','当','lambna=',num2str(lambna(k))))fprintf('近似解解析解误差')fprintf('\n')disp([u(:,N+1),w,ez])subplot(1,2,k)plot(x,u(:,N+1),'r*')holdon;plot(x,w)title(strcat('近似解与解析解比较',':','当','lambna=',num2str(lambna(k))))legend('近似解','解析解')end（2）functionkuosan2%向后差分X=0:0.05:1;U=exp(-pi^2*0.1)*sin(pi*X);x=0:0.1:1;[u1]=XH(0.1,100);[u2]=XH(0.5,20);plot(x,u1,'r:',x,u2,'g-.',X,U,'b','LineWidth',2);xlabel('x');ylabel('u');legend('\lambda=0.1','\lambda=0.5','准确解');title('向后差分');function[u]=XH(A,T)x=0:0.1:1;m=length(x);u=sin(pi*x);C1=sparse(1:m-2,1:m-2,1+2*A,m-2,m-2);C2=sparse(2:m-2,1:m-3,-A,m-2,m-2);C3=sparse(1:m-3,2:m-2,-A,m-2,m-2);C=C1+C2+C3;%系数矩阵forn=1:Tb=u(2:m-1);u(2:m-1)=C\b';end（3）functionkuosan3%C-N格式X=0:0.05:1;U=exp(-pi^2*0.1)*sin(pi*X);x=0:0.1:1;[u1]=CN(0.1,100);[u2]=CN(0.5,20);plot(x,u1,'r:',x,u2,'g-.',X,U,'b','LineWidth',2);xlabel('x');ylabel('u');legend('\lambda=0.1','\lambda=0.5','准确解');title('C-N格式');function[u]=CN(A,T)x=0:0.1:1;m=length(x);u=sin(pi*x);C1=sparse(1:m-2,1:m-2,1+A,m-2,m-2);C2=sparse(2:m-2,1:m-3,-0.5*A,m-2,m-2);C3=sparse(1:m-3,2:m-2,-0.5*A,m-2,m-2);C=C1+C2+C3;%系数矩阵forn=1:Tb=(1-A)*u(2:m-1)+(0.5*A)*u(3:m)+(0.5*A)*u(1:m-2);u(2:m-1)=C\b';end（4）functionkuosan4%DuFort-Frankel格式X=0:0.05:1;U=exp(-pi^2*0.1)*sin(pi*X);x=0:0.1:1;[u1]=DF(0.1,100);[u2]=DF(0.5,20);plot(x,u1,'r:',x,u2,'g-.',X,U,'b','LineWidth',2);xlabel('x');ylabel('u');legend('\lambda=0.1','\lambda=0.5','准确解');title('DuFort-Frankel格式');function[u]=DF(A,T)x=0:0.1:1;m=length(x);u=sin(pi*x);b=u;v=u;forn=1:Tfori=2:m-1u(i)=((1-2*A)*b(i)+2*A*v(i+1)+2*A*v(i-1))/(1+2*A);endb=v;v=u;end结论：显示格式向前格式和DuFort-Frankel格式的误差较小，向后差分格式的误差最大。