微分方程数值解第四次上机报告.doc

南京信息工程大学 实验（实习）报告 实验课程 微分方程数值解 实验名称 第四次实验 实验日期 2016 指导老师 专业 年级 姓名 学号 得分 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 实验目的： *编程五点三层格式求解两道二阶波动方程定解问题； *讨论不同网格划分得到的数值求解结果. 实验内容： 编程内容 对二阶波动方程定解问题（1）求解的程序见wave.m文件，对定解问题（2）求解的程序见wave2.m文件. 问题求解 *先求解如下一道二阶波动方程定解问题（t范围取有限）： （1） 数值求解： 用五点三层差分格式求解，先取空间步长,时间步长,故柯朗数,满足稳定性条件.输出数值解三维图和末时刻解如图1. 图1：疏网格时数值解三维图（左）和末时刻解图像（右） 从图1可以看出数值解图像纹理粗糙，末时刻解还有微量波动，这时可以加密网格获得更好的解.为加密网格，再取空间步长,时间步长.柯朗数不变，仍满足稳定性条件.输出数值解三维图和末时刻解如图2. 从图2可以看出数值解三维图纹理相对图1更为平滑，且直观上看末时刻解比较稳定。 图2：较密网格时数值解三维图（左）和末时刻解图像（右） 再改动网格划分，取空间步长,时间步长,柯朗数 ，仍满足稳定性条件. 输出数值解三维图和末时刻解如图3。 图3：较密网格时数值解三维图（左）和末时刻解图像（右） 从图3可以看出数值解三维图直观上看与图2几乎无差别，仍能很好刻画出精确解。且直观上看末时刻解依旧比较稳定。 继续改动网格划分，取空间步长,时间步长,柯朗数 ，此时将不满足稳定性条件. 输出数值解三维图和末时刻解如图4。 图4：较疏网格时数值解三维图（左）和末时刻解图像（右） 从图4可以看出当柯朗数不满足稳定性条件时，数值解三维图刻画效果很差。且末时刻解波动极大，很不稳定。 *再求解如下一道二阶波动方程定解问题（t范围取有限）： （2） 分离变量法易求得其解析解为： 数值求解： 用五点三层格式求解，具体程序见wave2.m文件。 精确解取60项，编程构造精确解矩阵： u_exactc=zeros(60,M,N); u_ec=zeros(M,N); for i=1:M for j=1:N for k=1:60 %取60项 u_exactc(k,i,j)=4/pi^4*(1-(-1)^(k))/(k^4)*sin(k*pi*t(j))*sin(k*pi*x(i)); u_ec(i,j)=sum(u_exactc(:,i,j)); end end end u_ec=u_ec; 再构造末时刻精确解： 同样是取60项，故可由上面构造的精确解矩阵最后行得出： u_e=u_ec(end,:); 考虑计算精确解有相当的时间复杂度，故这里使用较疏网格划分。取空间步长,时间步长，柯朗数满足稳定性条件。编程输出精确解三维图，数值解三维图，绝对误差三维图和末时刻的精确解及数值解图像如下： 图5（a）：精确解三维图 图5（b）：数值解三维图 图5（c）：绝对误差三维图 图5（d）：末时刻精确解与数值解 图5：问题（2）求解输出图像 从数值解和精确解图像上看，数值方法很好刻画出了精确解。从绝对误差上看，整个区域内的误差都在数量级，也说明方法比较稳定。从末时刻精确解与数值解图像上看，两者较为贴合，无法直观上看出差别，也说明了方法的稳定性。 *注：这里对范围进行限制，取，右端不取1。原因是末时刻精确解在处为零向量，而数值解就不是零向量，绘制出的末时刻精确解数值解图像两者不吻合（数值解与精确解整体三维图直观上仍然是吻合的）。如下图： 图6：取时末时刻精确解与解析解 容易理解图中数值解曲线即为绝对误差曲线。造成这种现象的原因实际上是因为坐标刻度不同。从上图可见图像纵坐标给出的范围在数量级，实则满足误差精度要求，但是图形不如图5（d）美观。故这里取。 程序附录 wave.m文件： clc clear %二阶波动方程求解 c=1; %系数c的值 a=1; %x范围 b=1; %t范围 dx=0.02; %x方向步长deltax dt=0.01