微分方程数值解第二次上机报告

1、实验(实习)报告实验日期2016指导老师学号得分南京信息工程大学实验课程微分方程数值解实验名称第二次实验专业年级姓名实验目的:* 对一维抛物方程(热传导方程)数值求解的向前差分格式、向后差分格式、Crank-Nicolson格式、交替显隐式格式以及跳点格式进行编程。* 求解一道抛物方程定解问题实例，编程给出解三维图以及末时刻精确解与数值解比较图。* 以理论推导和编程实现分析差分方法分别对空间步长及时间步长的整体误差精度(收敛阶)(以向前差分为例)。实验内容：1 .方法编程已编程向前差分格式、向后差分格式、CrankNicolson格式、交替显隐式格式以及跳点格式。分别见附录的forward.m

2、文件、backward.m文件、Crank_Nicolson.m文件、AFB.m文件和skip.m文件。以向前差分为例探讨差分方法分别对空间步长及时间步长的收敛阶的程序分别见附录中errforwardx.m文件和errforwardt.m文件。2 .实例求解下求解一个一维热传导方程定解问题实例：dud2u;1J：dtdx22x,0<x«u(x,0)=<0.1；口1-x),u(0,t)=u1,t)=00<x<1,0<t<用偏微分方程相关知识可以求得其解析解为：匕81sin2_2u(x,t)=一%(R'k)tsin(k二x).二2k=1k2解析

3、解参与编程时，其中k取60,可以达到精确解精度要求。分别运行五个差分格式，得出解三维图和末时刻精确解与格式数值解对比图由于五个格式运行结果无法从肉眼观察出区别，下取向前差分格式运行结果展示.如下图图1、图2和图3。向前差分求解MM又图1:向前差分求解结果三维图卡时制日前整分力翻:蓝琳解也般图图2:末时刻数值解与精确解比较图3:数值解与精确解比较放大图结果分析：从图1可以看出u(x,t)关于时间t逐渐降低，在边界x=0,x=1上满足零边界值条件，符合热传导物理规律。从图2可以看出向前差分格式在空间步长h=Ax=0.02，时间步长k=&=0.0002时，有很好的收敛性。末时刻差分方法解与精

4、确解较为贴合。经k1过多次放大后(如图3)才可以看出细微差别。这是因为此时网格比s=a2。=-,h22满足稳定性条件。而五种方法原理不同，取不同空间步长、时间步长时的收敛性也不同。下面研究差分方法的收敛阶，以深入研究收敛性。3 .收敛阶理论推导将方程U=auxx在节点(Xj,tn)离散化：-n2n孚1=a与Ij=1,2，J;n=1,2,，N.(1)口"x时间步长为k,空间步长为h.*下以向前差分差分格式为例：对充分光滑的解u,由Taylor展式：u(Xj,tn)k2?2u(Xj,tn)3、U(Xj,tn)=UXj,tn)+k；十+Qk).(2)ft2jtu(Xj.1,tn)=u(xj

5、,tn)-：U(Xj,tn)h3/u(Xj,tn)丫.x二u(Xj,tn)h2万2二U(Xj,tn):x2(3)u(Xj,tn)=u(Xj,tn)-h:x4:UXj,tn):XOh5).h22“Xj,tn):x2,23h:u(Xj,tn)3-3Xh4?u(Xj,tn)疝一千一Oh5).对(2)式移项得：L(Xj,tn)_u(Xj,tn1)-u(Xj,tn)_k："(Xj,tn)ft-k2ft2(3)、(4)式相加得：qk2).d2u(Xj,tn)u(Xj+,tn)-2u(Xj,tn)+u(Xj3)h2'(Xj,tn).22xh(5)、(6)式带入(1)式得：12fx4(4)(

6、5)qh3).(6)",tn1)7Xj,tn)二ak其中：u(Xj1,tn)-2u(Xj,tn)u(Xj.1,tn)n小、JRj(u).hRn(u).2k二u(Xj,tn)2ft2Ok2)qh3).(7)式舍去Rn(u)即得到逼近(1)式的向前格式差分方程:uf1-ujnu：1-2u;u：h,j=1,2,J;n=1,2,N.(8)(9)其中，ujn=u(Xj,tn).从(8)式来看，网格化近似后余项Rn(u)对时间步长k的局部误差阶为2,对空间步长h的局部误差阶为3.于是对时间、空间两种步长的整体误差阶为1和2.4 .收敛阶编程探讨以向前差分格式为例，先固定时间步长，变动空间步长：固

7、定时间步长k=0.00,分别取四个空间步长hi=0.0124h2=0.024h3=0.05,儿=0.1,这样取值是为了避免在右边界处数值计算时有时为t=l,有时取不到t=l,以影响末时刻结果精确性。计算末时刻相对误差，误差与步长分别取对数后绘图如图4。末时刻向前差分方法相对误差随空间中氏变化对数-行数图log(h)图4:末时刻向前差分方法相对误差随空间步长变化对数-对视图图中其实有三条直线，程序中用矩阵U存放了三条直线的斜率.分别约为：1.403、2.135、2.056.符合理论讨论中的关于空间步长达到2阶收敛精度。再固定空间步长，变动时间步长：固定空间步长h=0.01,取两个空间步长k1=0

8、.00001,k2=0.000011,再计算末时刻相对误差。误差与步长分别取对数后绘图如图5。同时计算了这条直线的斜率约为1.1593符合理论讨论中的关于时间步长达到1阶收敛精度。末时刻向刖差分方法相对误差随时间步长变化对数,对数图-784-786-7SS5O792794-7.96-11.52-11.5-11.40-11.46-11.44-1142-11.4log(k)图5:末时刻向前差分方法相对误差随时间步长变化对数-对视图5.附录所有Matlab程序如下:forward.m文件：clcclearl=1;%参数l的值a=1;%系数a的值tmax=0.1;%t最大值k=0.0002;%时间t增

9、量h=0.02;%x增量s=aA(2)*k/(hA2);%网格比x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);X,T=meshgrid(x,t);%u(x,0)初始层fori=1:oifx(i)<=1/2u(i,1)=2*x(i);elseu(i,1)=2-2*x(i);endend%u(0,t)和u(l,t)边界条件u(1,：)=0；u(o,：)=0；%向前差分forj=1:(p-1)fori=2:(o-1)u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endend%末时刻精

10、确解u_exact=zeros(60,o);fori=1:ofork=1:60%取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)A2)*sin(k*pi*x(i)/(piA2);u_e(i)=sum(u_exact(:,i);endend%解三维网格绘图figure()mesh(X',T',u)xlabel('x');ylabel('t');zlabel('u(x,t)');title('向前差分求解三维图')%末时刻解比较figure()plot(x,u(:,p),&#

11、39;-')%差分方法求解holdonplot(x,u_e,'-.')%精确解xlabel('x');ylabel('u(x,t)');title('末时刻向前差分方法解与精确解比较图')legend('向前差分方法解,'精确解(n=60)')backward.m文件：clcclearl=1;%参数l的值a=1;%系数a的值tmax=0.1;%t最大值k=0.00004;%时间t增量h=0.01;%x增量x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zer

12、os(o,p);s=aA(2)*k/(hA2);%网格比N=o-2;%每一层内点个数m,n=meshgrid(x,t);%u(x,0)初始层fori=1:oifx(i)<=1/2u(i,1)=2*x(i);elseu(i,1)=2-2*x(i);endend%u(0,t)和u(l,t)边界条件u(1,:)=0;u(o,:)=0;E=zeros(N,1);%边界不为零在此改动F=zeros(N,1);%隐式方程组右端先声明%隐式差分方程组系数矩阵AA=zeros(N,N);A(1,1)=1+2*s;A(1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;fori=2:N-1A

13、(i,i)=1+2*s;A(i,i+1)=-s;end%对时间层逐层求解forj=1:p-1F=u(2:N+1,j)+E;%方程组右端更新u(2:N+1,j+1)=AF;end%末时刻精确解u_exact=zeros(60,o);fori=1:ofork=1:60%取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)A2)*sin(k*pi*x(i)/(piA2);u_e(i)=sum(u_exact(:,i);endend%解三维网格绘图figure()mesh(m',n',u)xlabel('x');ylabel(&#

14、39;t');zlabel('u(x,t)');title('隐式差分求解三维图,)%末时刻解比较figure()plot(x,u(:,p),'-')%隐式差分方法求解holdonplot(x,u_e,'-.')%精确解xlabel('x');ylabel('u(x,t)');title('末时刻隐式差分方法解与精确解比较图')legend。隐式差分方法解,'精确解(n=60)')Crank_Nicolson.m文件：clcclearl=1;%参数l的值a=1;%系数

15、a的值tmax=0.1;%t最大值k=0.0002;%时间t增量h=0.02;%x增量x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);s=aA(2)*k/(2*hA2);%C-N网格比N=o-2;%每一层内点个数m,n=meshgrid(x,t);%u(x,0)初始层fori=1:oifx(i)<=1/2u(i,1)=2*x(i);elseu(i,1)=2-2*x(i);endend%u(0,t)和u(l,t)边界条件u(1,:)=0;u(o,:)=0;E=zeros(N,1);%非齐次边界问题在此改动F=zeros(N,1

16、);%隐式方程组右端先声明%方程组隐式一端系数矩阵AA=zeros(N,N);A(1,1)=1+2*s;A(1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;fori=2:N-1A(i,i)=1+2*s;A(i,i+1)=-s;A(i,i-1)=-s;end%方程组显式一端系数矩阵TT=zeros(N,N);T(1,1)=1-2*s;T(1,2)=s;T(N,N)=1-2*s;T(N,N-1)=s;fori=2:N-1T(i,i)=1-2*s;T(i,i+1)=s;T(i,i-1)=s;end%对时间层逐层求解forj=1:p-1F=u(2:N+1,j)+E;%方程组右端更新u

17、(2:N+1,j+1)=inv(A)*T*F;end%末时刻精确解u_exact=zeros(60,o);fori=1:ofork=1:60%取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)A2)*sin(k*pi*x(i)/(piA2);u_e(i)=sum(u_exact(:,i);endend%解三维网格绘图figure()mesh(m',n',u)xlabel('x');ylabel('t');zlabel('u(x,t)');title('C-N差分求解三维图,)%末时

18、刻解比较figure()plot(x,u(:,p),'-')%隐式差分方法求解holdonplot(x,u_e,'-.')%精确解xlabel('x');ylabel('u(x,t)');title(,末时刻C-N差分方法解与精确解比较图,)legend('C-N差分方法解,精确解(n=60),)AFB.m文件：clcclear%一维交替显隐式格式l=1;%参数l的值a=1;%系数a的值tmax=0.1;%t最大值k=0.0002;%时间t增量h=0.02;%x增量x=0:h:l;t=0:k:tmax;o=length(x

19、);p=length(t);u=zeros(o,2*p-1);%共2p-1层前p层为初始层+隐式矫正层后p-1层为显式预测层s=aA(2)*k/(2*hA2);%交替显隐式格式网格比N=o-2;%每一层内点个数m,n=meshgrid(x,t);%u(x,0)初始层fori=1:oifx(i)<=1/2u(i,1)=2*x(i);elseu(i,1)=2-2*x(i);endend%u(0,t)和u(l,t)边界条件u(1,:)=0;u(o,:)=0;E=zeros(N,1);%非齐次边界问题在此改动F=zeros(N,1);%隐式方程组右端先声明%方程组隐式一端系数矩阵AA=zeros

20、(N,N);A(1,1)=1+2*s;A(1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;fori=2:N-1A(i,i)=1+2*s;A(i,i+1)=-s;end%时间层循环forj=1:p-1fori=2:(o-1)u(i,j+p)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);%显式预估endF=u(2:N+1,j+p)+E;%方程组右端更新u(2:N+1,j+1)=AF;%隐式矫正end%末时刻精确解u_exact=zeros(60,o);fori=1:ofork=1:60%取60项u_exact(k,i)=8*sin(k*pi/2)*e

21、xp(-0.1*(k*pi)A2)*sin(k*pi*x(i)/(piA2);u_e(i)=sum(u_exact(:,i);endend%解三维网格绘图figure()mesh(m',n',u(:,1:p)xlabel('x');ylabel('t');zlabel('u(x,t)');title('交替显隐式格式求解三维图,)%末时刻解比较figure()plot(x,u(:,p),'-')%隐式差分方法求解holdonplot(x,u_e,'-.')%精确解xlabel('x&

22、#39;);ylabel('u(x,t)');title('末时刻交替显隐式格式解与精确解比较图')legend(交替显隐式格式解,'精确解(n=60)')skip.m文件:clcclear%一维跳点格式l=1;%参数l的值a=1;%系数a的值tmax=0.1;%t最大值k=0.0002;%时间t增量h=0.02;%x增量x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);s=aA(2)*k/(hA2);%跳点格式网格比N=o-2;%每一层内点个数m,n=meshgrid(x,t);%u

23、(x,0)初始层fori=1:oifx(i)<=1/2u(i,1)=2*x(i);elseu(i,1)=2-2*x(i);endend%u(0,t)和u(l,t)边界条件u(1,:)=0;u(o,:)=0;%跳点格式循环forj=1:p-1fori=2:(o-1)%(i+j+1)先在偶数格点显式%(i+j+1)再在奇数格点ifrem(i+j+1,2)=0u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endifrem(i+j+1,2)=0u(i,j+1)=u(i+1,j)/(1+2*s)+s*u(i+1,j+1)+s*u(i-1,j+1);隐式

24、(实际上显式)endendend%末时刻精确解u_exact=zeros(60,o);fori=1:ofork=1:60%取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)A2)*sin(k*pi*x(i)/(piA2);u_e(i)=sum(u_exact(:,i);endend%解三维网格绘图figure()mesh(m',n',u(:,1:p)xlabel('x');ylabel('t');zlabel('u(x,t)');title('跳点格式求解三维图,)%末时刻解比较f

25、igure()plot(x,u(:,p),'-')%跳点格式求解holdonplot(x,u_e,'-.')%精确解xlabel('x');ylabel('u(x,t)');title('末时刻跳点格式解与精确解比较图')legend(跳点格式解,'精确解(n=60)')errforwardx.m文件：clcclearformatlong%向前差分格式收敛阶l=1;%参数l的值a=1;%系数a的值tmax=0.1;%t最大值k=0.00004;%时间t增量h0=0.0125;%初始空间步长err=z

26、eros(1,4);%#放相对误差ke=1;m=log(1248*h0);forhh=1248%文寸空间步长循环h=hh*h0;%x增量s=aA(2)*k/(hA2);%网格比x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);%u(x,0)初始层fori=1:oifx(i)<=1/2u(i,1)=2*x(i);elseu(i,1)=2-2*x(i);endend%u(0,t)和u(l,t)边界条件u(1,:)=0;u(o,:)=0;%向前差分forj=1:(p-1)fori=2:(o-1)u(i,j+1)=s*u(i+1,j)

27、+(1-2*s)*u(i,j)+s*u(i-1,j);endend%末时刻精确解u_exact=zeros(60,o);u_e=zeros(1,o);forii=1:oforkk=1:60%取60项u_exact(kk,ii)=8*sin(kk*pi/2)*exp(-0.1*(kk*pi)A2)*sin(kk*pi*x(ii)/(piA2);u_e(ii)=sum(u_exact(:,ii);endenderr(ke)=norm(u(:,p)-u_e',2)/norm(u_e,2);%末时刻相对误差ke=ke+1;end;%各段斜率计算len=length(err);U=zeros(1,len-1);forflag=1:len-1U(flag)=(log(err(flag+1)-log(err(flag)/(m(flag+1)-m(flag);end%末时刻相对误差plot(m,log(err),'-')xlabel('log(h)');ylabel('log(error)');title('末时刻向前差分方法相对误差随空间步长变化对数-对数图')legend('向