资源描述：
《微分方程数值解第二次上机报告.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、南京信息工程大学实验（实习）报告实验课程微分方程数值解实验名称第二次实验实验日期2016指导老师专业年级姓名学号得分------------------------------------------------------------------实验目的：*对一维抛物方程（热传导方程）数值求解的向前差分格式、向后差分格式、Crank—Nicolson格式、交替显隐式格式以及跳点格式进行编程。*求解一道抛物方程定解问题实例,编程给出解三维图以及末时刻精确解与数值解比较图。*以理论推导和编程实现分析差分方法分别对空间步长及时间步长的
2、整体误差精度（收敛阶）（以向前差分为例）。实验内容：1.方法编程已编程向前差分格式、向后差分格式、Crank—Nicolson格式、交替显隐式格式以及跳点格式。分别见附录的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、时网格比,满足稳定性条件。而五种方法原理不同，取不同空间步长、时间步长时的收敛性也不同。下面研究差分方法的收敛阶，以深入研究收敛性。1.收敛阶理论推导将方程在节点离散化：（1）时间步长为，空间步长为.*下以向前差分差分格式为例：对充分光滑的解，由Taylor展式：（2）（3）（4）对（2）式移项得：（5）（3）、（4）式相加得：（6）（5）、（6）式带入（1）式得：（7）其中：（8）（7）式舍去即得到逼近（1）式的向前格式差分方程：（9）其中，从（8）式来看，网格化近似后余项对时间步长的局部误差阶为2，对空间步长的局部误差阶为3.
5、于是对时间、空间两种步长的整体误差阶为1和2.1.收敛阶编程探讨以向前差分格式为例，先固定时间步长，变动空间步长：固定时间步长，分别取四个空间步长，这样取值是为了避免在右边界处数值计算时有时为，有时取不到，以影响末时刻结果精确性。计算末时刻相对误差，误差与步长分别取对数后绘图如图4。图4：末时刻向前差分方法相对误差随空间步长变化对数-对视图图中其实有三条直线，程序中用矩阵U存放了三条直线的斜率.分别约为：1.403、2.135、2.056.符合理论讨论中的关于空间步长达到2阶收敛精度。再固定空间步长，变动时间步长：固定空间步长，取
6、两个空间步长，再计算末时刻相对误差。误差与步长分别取对数后绘图如图5。同时计算了这条直线的斜率约为符合理论讨论中的关于时间步长达到1阶收敛精度。图5：末时刻向前差分方法相对误差随时间步长变化对数-对视图1.附录所有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
7、(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
8、(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)xla
