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

1、南京信息工程大学 实验（实习）报告实验课程 微分方程数值解 实验名称 第二次实验 实验日期 2016 指导老师 专业 年级 姓名 学号 得分 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 实验目的：*对一维抛物方程（热传导方程）数值求解的向前差分格式、向后差分格式、CrankNicolson格式、交替显隐式格式以及跳点格式进行编程。*求解一道抛物方程定解问题实例,编程给出解三维图

2、以及末时刻精确解与数值解比较图。*以理论推导和编程实现分析差分方法分别对空间步长及时间步长的整体误差精度（收敛阶）（以向前差分为例）。实验内容：1. 方法编程已编程向前差分格式、向后差分格式、CrankNicolson格式、交替显隐式格式以及跳点格式。分别见附录的forward.m文件、backward.m文件、Crank_Nicolson.m文件、AFB.m文件和skip.m文件。以向前差分为例探讨差分方法分别对空间步长及时间步长的收敛阶的程序分别见附录中errforwardx.m文件和errforwardt.m文件。2. 实例求解下求解一个一维热传导方程定解问题实例：用偏微分方程相关知识可

3、以求得其解析解为：解析解参与编程时，其中k取60,可以达到精确解精度要求。分别运行五个差分格式，得出解三维图和末时刻精确解与格式数值解对比图。由于五个格式运行结果无法从肉眼观察出区别，下取向前差分格式运行结果展示。如下图图1、图2和图3。 图1：向前差分求解结果三维图图2：末时刻数值解与精确解比较图3：数值解与精确解比较放大图结果分析：从图1可以看出关于时间逐渐降低，在边界上满足零边界值条件，符合热传导物理规律。从图2可以看出向前差分格式在空间步长时间步长时，有很好的收敛性。末时刻差分方法解与精确解较为贴合。经过多次放大后（如图3）才可以看出细微差别。这是因为此时网格比,满足稳定性条件。而五种

4、方法原理不同，取不同空间步长、时间步长时的收敛性也不同。下面研究差分方法的收敛阶，以深入研究收敛性。3. 收敛阶理论推导 将方程在节点离散化： （1）时间步长为，空间步长为.*下以向前差分差分格式为例：对充分光滑的解，由Taylor展式：（2） （3） （4）对（2）式移项得： （5）（3）、（4）式相加得：（6）（5）、（6）式带入（1）式得： （7）其中： （8）（7）式舍去即得到逼近（1）式的向前格式差分方程： （9）其中，从（8）式来看，网格化近似后余项对时间步长的局部误差阶为2，对空间步长的局部误差阶为3.于是对时间、空间两种步长的整体误差阶为1和2.4. 收敛阶编程探讨以向前差分格

5、式为例，先固定时间步长，变动空间步长：固定时间步长，分别取四个空间步长，这样取值是为了避免在右边界处数值计算时有时为，有时取不到，以影响末时刻结果精确性。计算末时刻相对误差，误差与步长分别取对数后绘图如图4。图4：末时刻向前差分方法相对误差随空间步长变化对数-对视图图中其实有三条直线，程序中用矩阵U存放了三条直线的斜率.分别约为：1.403、2.135、2.056.符合理论讨论中的关于空间步长达到2阶收敛精度。再固定空间步长，变动时间步长：固定空间步长，取两个空间步长，再计算末时刻相对误差。误差与步长分别取对数后绘图如图5。同时计算了这条直线的斜率约为符合理论讨论中的关于时间步长达到1阶收敛精

6、度。图5：末时刻向前差分方法相对误差随时间步长变化对数-对视图5. 附录所有Matlab程序如下：forward.m文件：clcclearl=1; %参数l的值a=1; %系数a的值tmax=0.1; %t最大值k=0.0002; %时间t增量h=0.02; % x增量s=a(2)*k/(h2); %网格比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) 初始层for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(i,1)=2-2*x(i); e

7、ndend% u(0,t)和u(l,t) 边界条件u(1,:)=0;u(o,:)=0;% 向前差分for j=1:(p-1) for i=2:(o-1) u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endend% 末时刻精确解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)2)*sin(k*pi*x(i)/(pi2);u_e(i)=sum(u_exact(:,i); endend% 解三维网格绘图figure()mes

8、h(X,T,u)xlabel(x);ylabel(t);zlabel(u(x,t);title(向前差分求解三维图)% 末时刻解比较figure()plot(x,u(:,p),-) %差分方法求解hold onplot(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:t

9、max;o=length(x);p=length(t);u=zeros(o,p);s=a(2)*k/(h2); %网格比N=o-2; %每一层内点个数m,n=meshgrid(x,t);% u(x,0) 初始层for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(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(

10、1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;for i=2:N-1 A(i,i)=1+2*s; A(i,i+1)=-s; A(i,i-1)=-s;end% 对时间层逐层求解for j=1:p-1 F=u(2:N+1,j)+E; %方程组右端更新 u(2:N+1,j+1)=AF;end% 末时刻精确解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)2)*sin(k*pi*x(i)/(pi2);u_e(i)=sum(u_exact(:,i); end

11、end% 解三维网格绘图figure()mesh(m,n,u)xlabel(x);ylabel(t);zlabel(u(x,t);title(隐式差分求解三维图)% 末时刻解比较figure()plot(x,u(:,p),-) %隐式差分方法求解hold onplot(x,u_e,-.) %精确解xlabel(x);ylabel(u(x,t);title(末时刻隐式差分方法解与精确解比较图)legend(隐式差分方法解,精确解（n=60）)Crank_Nicolson.m文件：clcclearl=1; %参数l的值a=1; %系数a的值tmax=0.1; %t最大值k=0.0002; %时间t

12、增量h=0.02; % x增量x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);s=a(2)*k/(2*h2); %C-N网格比N=o-2; %每一层内点个数m,n=meshgrid(x,t);% u(x,0) 初始层for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(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); %隐式方程组右端 先声明%方程组

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

14、 u(2:N+1,j+1)=inv(A)*T*F;end% 末时刻精确解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)2)*sin(k*pi*x(i)/(pi2);u_e(i)=sum(u_exact(:,i); endend% 解三维网格绘图figure()mesh(m,n,u)xlabel(x);ylabel(t);zlabel(u(x,t);title(C-N差分求解三维图)% 末时刻解比较figure()plot(x,u(:,p),-) %隐式差分方法求解h

15、old onplot(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);p=length(t);u=zeros(o,2*p-1); %共2p-1层 前p层为初始层+隐式矫正层 后p-1层为显式预测层s=a(2)*k/(2*h2

16、); %交替显隐式格式网格比N=o-2; %每一层内点个数m,n=meshgrid(x,t);% u(x,0) 初始层for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(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); %隐式方程组右端 先声明%方程组 隐式一端系数矩阵A A=zeros(N,N); A(1,1)=1+2*s;A(1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;for i=2:

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

18、)/(pi2);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),-) %隐式差分方法求解hold onplot(x,u_e,-.) %精确解xlabel(x);ylabel(u(x,t);title(末时刻交替显隐式格式解与精确解比较图)legend(交替显隐式格式解,精确解（n=60）)skip.m文件：clcclear%一维跳点格式l=1; %参

19、数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=a(2)*k/(h2); %跳点格式网格比N=o-2; %每一层内点个数m,n=meshgrid(x,t);% u(x,0) 初始层for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(i,1)=2-2*x(i); endend% u(0,t)和u(l,t) 边界条件u(1,:)=0;u(o,:)=0;%跳点格式循环for j=1

20、:p-1 for i=2:(o-1) if rem(i+j+1,2)=0 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); %(i+j+1)先在偶数格点显式 end if rem(i+j+1,2)=0 u(i,j+1)=u(i+1,j)/(1+2*s)+s*u(i+1,j+1)+s*u(i-1,j+1); %(i+j+1)再在奇数格点隐式（实际上显式） end end end% 末时刻精确解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60项u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1

21、*(k*pi)2)*sin(k*pi*x(i)/(pi2);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),-) %跳点格式求解hold onplot(x,u_e,-.) %精确解xlabel(x);ylabel(u(x,t);title(末时刻跳点格式解与精确解比较图)legend(跳点格式解,精确解（n=60）)errforwardx.m文件:clcc

22、learformat long%向前差分格式收敛阶l=1; %参数l的值a=1; %系数a的值tmax=0.1; %t最大值k=0.00004; %时间t增量h0=0.0125; %初始空间步长err=zeros(1,4); %存放相对误差ke=1;m=log( 1 2 4 8*h0); for hh= 1 2 4 8 %対空间步长循环 h=hh*h0; % x增量s=a(2)*k/(h2); %网格比x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);% u(x,0) 初始层for i=1:o if x(i)=1/2 u(i,1)=

23、2*x(i); else u(i,1)=2-2*x(i); endend% u(0,t)和u(l,t) 边界条件u(1,:)=0;u(o,:)=0;% 向前差分for j=1:(p-1) for i=2:(o-1) u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endend% 末时刻精确解u_exact=zeros(60,o);u_e=zeros(1,o);for ii=1:o for kk=1:60 %取60项u_exact(kk,ii)=8*sin(kk*pi/2)*exp(-0.1*(kk*pi)2)*sin(kk*pi*x(ii)/(pi

24、2);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);for flag=1:len-1 U(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(向前差分方法)errforwardt.m文件:clcclearformat longl=1; %参数l的值a=1; %系数a的值tmax=0.1; %t最大值h=0.01; %x增量k0=0.00