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

南京信息工程大学 实验（实习）报告 实验课程 微分方程数值解 实验名称 第二次实验 实验日期 2016 指导老师 专业 年级 姓名 学号 得分 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 实验目的： *对一维抛物方程（热传导方程）数值求解的向前差分格式、向后差分格式、Crank—Nicolson格式、交替显隐式格式以及跳点格式进行编程。 *求解一道抛物方程定解问题实例,编程给出解三维图以及末时刻精确解与数值解比较图。 *以理论推导和编程实现分析差分方法分别对空间步长及时间步长的整体误差精度（收敛阶）（以向前差分为例）。 实验内容： 方法编程 已编程向前差分格式、向后差分格式、Crank—Nicolson格式、交替显隐式格式以及跳点格式。分别见附录的forward.m文件、backward.m文件、Crank_Nicolson.m文件、AFB.m文件和skip.m文件。 以向前差分为例探讨差分方法分别对空间步长及时间步长的收敛阶的程序分别见附录中errforwardx.m文件和errforwardt.m文件。 实例求解 下求解一个一维热传导方程定解问题实例： 用偏微分方程相关知识可以求得其解析解为： 解析解参与编程时，其中k取60,可以达到精确解精度要求。 分别运行五个差分格式，得出解三维图和末时刻精确解与格式数值解对比图。 由于五个格式运行结果无法从肉眼观察出区别，下取向前差分格式运行结果展示。如下图图1、图2和图3。 图1：向前差分求解结果三维图 图2：末时刻数值解与精确解比较 图3：数值解与精确解比较放大图 结果分析： 从图1可以看出关于时间逐渐降低，在边界上满足零边界值条件，符合热传导物理规律。 从图2可以看出向前差分格式在空间步长时间步长时，有很好的收敛性。末时刻差分方法解与精确解较为贴合。经过多次放大后（如图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。 图4：末时刻向前差分方法相对误差随空间步长变化对数-对视图 图中其实有三条直线，程序中用矩阵U存放了三条直线的斜率.分别约为：1.403、2.135、2.056.符合理论讨论中的关于空间步长达到2阶收敛精度。 再固定空间步长，变动时间步长： 固定空间步长，取两个空间步长，再计算末时刻相对误差。误差与步长分别取对数后绘图如图5。 同时计算了这条直线的斜率约为符合理论讨论中的关于时间步长达到1阶收敛精度。 图5：末时刻向前差分方法相对误差随时间步长变化对数-对视图 附录 所有Matlab程序如下： forward.m文件： clc clear l=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) 初始层 for i=1:o if x(i)