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

南京信息工程大学实验（实习）报告实验课程微分方程数值解实验名称第二次实验实验日期指导老师王曰朋专业数学与应用数学年级姓名学号得分------------------------------------------------------------------实验目的：*对一维抛物方程（热传导方程）数值求解的向前差分格式、向后差分格式、Crank—Nicolson格式、交替显隐式格式以及跳点格式进行编程。*求解一道抛物方程定解问题实例,编程给出解三维图以及末时刻精确解与数值解比较图。*以理论推导和编程实现分析差分方法分别对空间步长及时间步长的整体误差精度（以向前差分为例）。实验内容：1.方法编程编程向前差分格式、向后差分格式、Crank—Nicolson格式、交替显隐式格式以及跳点格式。分别见附录的forward.m文件、backward.m文件、Crank_Nicolson.m文件、AFB.m文件和skip.m文件。以向前差分为例探讨差分方法分别对空间步长及时间步长的收敛阶的程序分别见附录。2.实例求解求解一个一维热传导方程定解问题实例：;1.00,10;0),1(),0(;121),1(2;210,2)0,(;22txtutuxxxxxudxuddtdu用偏微分方程相关知识可以求得其解析解为：1)(22).sin()2sin(8),(2ktkxkekktxu解析解参与编程时，其中k取60,可以达到精确解精度要求。分别运行五个差分格式，得出解三维图和末时刻精确解与格式数值解对比图。由于五个格式运行结果无法从肉眼观察出区别，下取向前差分格式运行结果展示。如下图图1、图2和图3。图1：向前差分求解结果三维图图2：末时刻数值解与精确解比较图3：数值解与精确解比较放大图结果分析：从图1可以看出),(txu关于时间t逐渐降低，在边界1,0xx上满足零边界值条件，符合热传导物理规律。从图2可以看出向前差分格式在空间步长,02.0xh时间步长0002.0tk时，有很好的收敛性。末时刻差分方法解与精确解较为贴合。经过多次放大后（如图3）才可以看出细微差别。这是因为此时网格比NeumannVonhkas满足,2122稳定性条件。而五种方法原理不同，取不同空间步长、时间步长时的收敛性也不同。下面研究这五种方法的收敛阶，以深入研究收敛性。3.收敛阶理论推导将方程xxtauu在节点),(njtx离散化：.,,2,1;,,2,1;22NnJjxuatunjnj（1）时间步长为k，空间步长为h.*下以向前差分差分格式为例：对充分光滑的解u，由Taylor展式：).(),(2),(),(),(32221kOttxukttxuktxutxunjnjnjnj（2）).(),(!4),(!3),(2),(),(),(54443332221hOxtxuhxtxuhxtxuhxtxuhtxutxunjnjnjnjnjnj（3）).(),(!4),(!3),(2),(),(),(54443332221hOxtxuhxtxuhxtxuhxtxuhtxutxunjnjnjnjnjnj（4）对（2）式移项得：).(),(2),(),(),(2221kOttxukktxutxuttxunjnjnjnj（5）（3）、（4）式相加得：).(),(12),(),(2),(),(344221122hOxtxuhhtxutxutxuxtxunjnjnjnjnj（6）（5）、（6）式带入（1）式得：).(),(),(2),(),(),(2111uRhtxutxutxuaktxutxunjnjnjnjnjnj（7）其中：).()(),(12),(2)(3244222hOkOxtxuhttxukuRnjnjnj（8）（7）式舍去)(uRnj即得到逼近（1）式的向前格式差分方程：.,,2,1;,,2,1,22111NnJjhuuuakuunjnjnjnjnj（9）其中，).,(njnjtxuu从（8）式来看，网格化近似后余项)(uRnj对时间步长k的局部误差阶为2，对空间步长h的局部误差阶为3.于是对时间、空间两种步长的整体误差阶为1和2.4.收敛阶编程探讨以向前差分格式为例，先固定时间步长，变动空间步长：固定时间步长0.00004k，取两个空间步长023.0,022.021hh，再计算末时刻相对误差。误差与步长分别取对数后绘图如下。图4：末时刻向前差分方法相对误差随空间步长变化对数-对视图同时计算了这条直线的斜率约为.1475.2符合理论讨论中的关于空间步长达到2阶收敛精度。再固定空间步长，变动时间步长：固定空间步长0.01h，取两个空间步长0.000011k0.00001,21k，再计算末时刻相对误差。误差与步长分别取对数后绘图如下。图5：末时刻向前差分方法相对误差随时间步长变化对数-对视图同时计算了这条直线的斜率约为1593.1符合理论讨论中的关于时间步长达到1阶收敛精度。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/(h^2);%网格比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%末时刻精确解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)^2)*sin(k*pi*x(i))/(pi^2);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),'-')%差分方法求解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=zeros(o,p);s=a^(2)*k/(h^2);%网格比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(i,i)=1+2*s;A(i,i+1)=-s;A(i,i-1)=-s;end%对时间层逐层求解forj=1:p-1F=u(2:N+1,j)+E;%方程组右端更新u(2:N+1,j+1)=A\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)^2)*sin(k*pi*x(i))/(pi^2);u_e(i)=sum(u_exact(:,i));endend%解三维网格绘图figure()mesh(m',n',u)xlabel('x');ylabel('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;%系数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/(2*h^2);%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);%隐式方程组右端先声明%方程组隐式一端系数矩阵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(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)^2)*sin(k*pi*x(i))/(pi^2);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),'-')%隐式差分方法求解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