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

南京信息工程大学 实验（实习）报告
实验课程 微分方程数值解 实验名称 第四次实验 实验日期 指导老师 王曰朋 专业 数学与应用数学 年级 姓名 学号 得分
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
实验目的：
*编程五点三层格式求解两道二阶波动方程定解问题； *讨论不同网格划分得到的数值求解结果. 实验内容： 1. 编程内容
对二阶波动方程定解问题（1）求解的程序见wave.m 文件，对定解问题（2）求解的程序见wave2.m 文件.
2. 问题求解
*先求解如下一道二阶波动方程定解问题（t 范围取有限）：
⎪⎩
⎪
⎨⎧≤≤==≤≤==<<<<=--.10,0)0,(,)0,(,10,
0),1(),0(10,10,u 2)3.0(400x x u e x u t t u t u t x u t x xx tt （1） 数值求解：
用五点三层差分格式求解，先取空间步长02.0=∆x ,时间步长01.0=∆t ,故柯
朗数15.0<=∆∆=x
t
c r ,满足稳定性条件.输出数值解三维图和末时刻解如图
1.
图1：疏网格)5.0(=r 时数值解三维图（左）和末时刻解图像（右）
从图1可以看出数值解图像纹理粗糙，末时刻解还有微量波动，这时可以加密网格获得更好的解.为加密网格，再取空间步长002.0=∆x ,时间步长001.0=∆t .柯朗数不变，仍满足稳定性条件.输出数值解三维图和末时刻解如图2.
从图2可以看出数值解三维图纹理相对图1更为平滑，且直观上看末时刻解比较稳定。

图2：较密网格)5.0(=r 时数值解三维图（左）和末时刻解图像（右） 再改动网格划分，取空间步长001.0=∆x ,时间步长001.0=∆t ,柯朗数 11≤=∆∆=x
t c r ，仍满足稳定性条件. 输出数值解三维图和末时刻解如图3。

图3：较密网格)1(=r 时数值解三维图（左）和末时刻解图像（右）
从图3可以看出数值解三维图直观上看与图2几乎无差别，仍能很好刻画出精确解。

且直观上看末时刻解依旧比较稳定。

继续改动网格划分，取空间步长01.0=∆x ,时间步长012.0=∆t ,柯朗数 12.1>=∆∆=x
t c r ，此时将不满足稳定性条件. 输出数值解三维图和末时刻解如
图4。

图4：较疏网格)2.1(=r 时数值解三维图（左）和末时刻解图像（右） 从图4可以看出当柯朗数不满足稳定性条件时，数值解三维图刻画效果很差。

且末时刻解波动极大，很不稳定。

*再求解如下一道二阶波动方程定解问题（t 范围取有限）：
⎪⎩⎪
⎨⎧≤≤-==≤≤==<<<<=.10),1()0,(,0)0,(,1.10,
0),1.1(),0(1.10,10,u x x x x u x u t t u t u t x u t xx tt （2） 分离变量法易求得其解析解为：
.sin sin )1(14),(1
4
4x n t n n
t x u n n
πππ
∑∞=--= 数值求解：
用五点三层格式求解，具体程序见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,:);
考虑计算精确解有相当的时间复杂度，故这里使用较疏网格划分。

取空间步长02.0=∆x ,时间步长01.0=∆t ，柯朗数满足稳定性条件。

编程输出精确解三维图，数值解三维图，绝对误差三维图和末时刻的精确解及数值解图像如下：
图5（a ）：精确解三维图
图5（b ）：数值解三维图
图5（c ）：绝对误差三维图
图5（d ）：末时刻精确解与数值解
图5：问题（2）求解输出图像
从数值解和精确解图像上看，数值方法很好刻画出了精确解。

从绝对误差上看，整个区域内的误差都在510-数量级，也说明方法比较稳定。

从末时刻精确解与数值解图像上看，两者较为贴合，无法直观上看出差别，也说明了方法的稳定性。

*注：这里对t 范围进行限制，取1.10≤≤t ，右端不取1。

原因是末时刻精确解在1=t 处为零向量，而数值解就不是零向量，绘制出的末时刻精确解数值解图像两者不吻合（数值解与精确解整体三维图直观上仍然是吻合的）。

如下图：
图6：取1
≤t时末时刻精确解与解析解
0≤
容易理解图中数值解曲线即为绝对误差曲线。

造成这种现象的原因实际上是因为坐标刻度不同。

从上图可见图像纵坐标给出的范围在5
10-数量级，实则满足误差精度要求，但是图形不如图5（d）美观。

故这里取1.1
≤t。

0≤。