二维波动方程的有限差分法

学 生 实 验 报 告实验课程名称 偏微分方程数值解 开课实验室 数统学院 学 院 数 统 年级 2013 专业班 信计02班 学 生 姓 名 学 号 开 课 时 间 2015 至 2016学年第 2 学期总 成 绩教师签名数学与统计学院制开课学院、实验室： 数统学院 实验时间 ： 2016年 6月20日实验项目名 称二维波动方程的有限差分法实验项目类型验证演示综合设计其他指导教师曾芳成 绩是一实验目的通过该实验，要求学生掌握求解二维波动方程的有限差分法，并能通过计算机语言编程实现。二实验内容考虑如下的初值问题： (1)1在第三部分写出问题（1）三层显格式。2根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。3取，分别将时刻的数值解画图显示。4. 该问题的解析解为，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。三实验原理、方法（算法）、步骤网格划分，故，。在内网点，利用二阶中心差商，对（1）建立差分格式： （2）整理得到： （3）其中，网比，局部截断误差为。考虑边界条件，差分格式为： （4）考虑初始条件，差分格式为： （5）考虑初始条件，利用二阶差商近似： （6）设时刻的点为内点，则满足差分格式（2），代入上式得到： （7）将（6）得到的结果代入（7）中，整理得到： （8）综上（2）、（4）、（5）、（8）得到三层显格式的差分格式为： （9）其中，局部截断误差为。四实验环境（所用软件、硬件等）及实验数据文件Matlab%二维波动方程数值计算（关键：怎么运用i,j,k三个指标建立循环）clc;%可以将代码换成函数m文件h=0.1;tau=0.1*h;%定义步长r=tau/h;%网比x,y,t=meshgrid(0:h:1,0:h:1,0:tau:1.4);%空间网格剖分uu=cos(sqrt(2)*pi*t).*sin(pi*x).*sin(pi*y);%精确解计算%第一层网点计算u=sin(pi*x).*sin(pi*y);%初始条件u1=u(:,:,1);%因为此时得到的u为11x11x141，故只取第一层%第二层网点计算for i=2:10 for j=2:10 u(i,j,2)=0.5*r2*(u(i+1,j,1)+u(i-1,j,1)+u(i,j+1,1)+u(i,j-1,1)+(1-2*r2)*u(i,j,1); u(11,:,2)=0;u(:,11,2)=0; endendu2=u(:,:,2);%第3-141层网点计算for k=2:140 for i=2:10 for j=2:10 u(i,j,k+1)=r2*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k)+(2-4*r2)*u(i,j,k)-u(i,j,k-1); u(11,:,k+1)=0;u(:,11,k+1)=0; end endend%结果分析与作图%wucha=abs(u-uu);%求绝对误差矩阵11x11x141wucha1=wucha(:,:,11);%计算t=0.1时刻的绝对误差矩阵11x11wucha2=wucha(:,:,51);%计算t=0.5时刻的绝对误差矩阵11x11wucha3=wucha(:,:,101);%计算t=1.0时刻的绝对误差矩阵11x11wucha4=wucha(:,:,141);%计算t=1.4时刻的绝对误差矩阵11x11x0=0:h:1;y0=0:h:1;%误差分析%作t=0.1时刻的绝对误差图subplot(2,2,1);mesh(x0,y0,wucha1);title(t=0.1时刻的绝对误差);xlabel(x变量);ylabel(y变量);zlabel(绝对误差值);%作t=0.5时刻的绝对误差图subplot(2,2,2);mesh(x0,y0,wucha2);title(t=0.5时刻的绝对误差);xlabel(x变量);ylabel(y变量);zlabel(绝对误差值);%作t=1.0时刻的绝对误差图subplot(2,2,3);mesh(x0,y0,wucha3);title(t=1.0时刻的绝对误差);xlabel(x变量);ylabel(y变量);zlabel(绝对误差值);%作t=1.4时刻的绝对误差图subplot(2,2,4);mesh(x0,y0,wucha4);title(t=1.4时刻的绝对误差);xlabel(x变量);ylabel(y变量);zlabel(绝对误差值);%四个时刻数值解、精确解%作t=0.1、0.5时刻的数值解与精确解subplot(2,2,1);mesh(x0,y0,u(:,:,11);%作t=0.1时刻的数值解title(t=0.1时刻的数值解);xlabel(x变量);ylabel(y变量);zlabel(u值);subplot(2,2,2);mesh(x0,y0,uu(:,:,11);%作t=0.1时刻的精确解title(t=0.1时刻的精确解);xlabel(x变量);ylabel(y变量);zlabel(u值);%作t=0.5时刻的数值解与精确解subplot(2,2,3);mesh(x0,y0,u(:,:,51);%作t=0.5时刻的数值解title(t=0.5时刻的数值解);xlabel(x变量);ylabel(y变量);zlabel(u值);subplot(2,2,4);mesh(x0,y0,uu(:,:,51);%作t=0.5时刻的精确解title(t=0.5时刻的精确解);xlabel(x变量);ylabel(y变量);zlabel(u值);%分别复制粘贴运行%作t=1.0、1.4时刻的数值解与精确解subplot(2,2,1);mesh(x0,y0,u(:,:,101);%作t=1.0时刻的数值解title(t=1.0时刻的数值解);xlabel(x变量);ylabel(y变量);zlabel(u值);subplot(2,2,2);mesh(x0,y0,uu(:,:,101);%作t=1.0时刻的精确解title(t=1.0时刻的精确解);xlabel(x变量);ylabel(y变量);zlabel(u值);%作t=1.4时刻的数值解与精确解subplot(2,2,3);mesh(x0,y0,u(:,:,141);%作t=1.4时刻的数值解title(t=1.4时刻的数值解);xlabel(x变量);ylabel(y变量);zlabel(u值);subplot(2,2,4);mesh(x0,y0,uu(:,:,141);%作t=1.4时刻的精确解title(t=1.4时刻的精确解);xlabel(x变量);ylabel(y变量);zlabel(u值);五实验结果及实例分析1、时刻的数值解与精确解图 图1 t=0.1、0.5时刻的数值解、精确解 图2 t=1.0、1.4时刻的数值解、精确解注：上两图为四个时刻的数值解与精确解，三层显格式达二阶收敛，不难看出，收敛效果很好，符合理论。下图是四个时刻的绝对误差图像，从图中看出，绝对误差较小，且经过计算得到，收敛阶近似于2，正好符合理论值。2、时刻的绝对误差图图3 四个时刻的绝对误差3、四个时刻（t=0.1、0.5、1.0、1.4）的绝对误差表t=0.1时刻的绝对误差0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0002 0.0001 0.0001 0.0000 0.0000 0.0001 0.0003 0.0004 0.0004 0.0005 0.0004 0.0004 0.0003 0.0001 0.0000 0.0000 0.0002 0.0004 0.0005 0.0006 0.0006 0.0006 0.0005 0.0004 0.0002 0.0000 0.0000 0.0002 0.0004 0.0006 0.0007 0.0007 0.0007 0.0006 0.0004 0.0002 0.0000 0.0000 0.0002 0.0005 0.0006 0.0007 0.0008 0.0007 0.0006 0.0005 0.0002 0.0000 0.0000 0.0002 0.0004 0.0006 0.0007 0.0007 0.0007 0.0006 0.0004 0.0002 0.0000 0.0000 0.0002 0.0004 0.0005 0.0006 0.0006 0.0006 0.0005 0.0004 0.0002 0.0000 0.0000 0.0001 0.0003 0.0004 0.0004 0.0005 0.0004 0.0004 0.0003 0.0001 0.0000 0.0000 0.0001 0.0001 0.0002 0.0002 0.0002 0.0002 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 t=0.5时刻的绝对误差0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0007 0.0013 0.0018 0.0021 0.0022 0.0021 0.0018 0.0013 0.0007 0.0000 0.0000 0.0013 0.0025 0.0034 0.0040 0.0042 0.0040 0.0034 0.0025 0.0013 0.0000 0.0000 0.0018 0.0034 0.0047 0.0055 0.0058 0.0055 0.0047 0.0034 0.0018 0.0000 0.0000 0.0021 0.0040 0.0055 0.0065 0.0068 0.0065 0.0055 0.0040 0.0021 0.0000 0.0000 0.0022 0.0042 0.0058 0.0068 0.0071 0.0068 0.0058 0.0042 0.0022 0.0000 0.0000 0.0021 0.0040 0.0055 0.0065 0.0068 0.0065 0.0055 0.0040 0.0021 0.0000 0.0000 0.0018 0.0034 0.0047 0.0055 0.0058 0.0055 0.0047 0.0034 0.0018 0.0000 0.0000 0.0013 0.0025 0.0034 0.0040 0.0042 0.0040 0.0034 0.0025 0.0013 0.0000 0.0000 0.0007 0.0013 0.0018 0.0021 0.0022 0.0021 0.0018 0.0013 0.0007 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 t=1.0时刻的绝对误差0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0016 0.0031 0.0043 0.0051 0.0053 0.0051 0.0043 0.0031 0.0016 0.0000 0.0000 0.0031 0.0059 0.0082 0.0096 0.0101 0.0096 0.0082 0.0059 0.0031 0.0000 0.0000 0.0043 0.0082 0.0113 0.0132 0.0139 0.0132 0.0113 0.0082 0.0043 0.0000 0.0000 0.0051 0.0096 0.0132 0.0156 0.0164 0.0156 0.0132 0.0096 0.0051 0.0000 0.0000 0.0053 0.0101 0.0139 0.0164 0.0172 0.0164 0.0139 0.0101 0.0053 0.0000 0.0000 0.0051 0.0096 0.0132 0.0156 0.0164 0.0156 0.0132 0.0096 0.0051 0.0000 0.0000 0.0043 0.0082 0.0113 0.0132 0.0139 0.0132 0.0113 0.0082 0.0043 0.0000 0.0000 0.0031 0.0059 0.0082 0.0096 0.0101 0.0096 0.0082 0.0059 0.0031 0.0000 0.0000 0.0016 0.0031 0.0043 0.0051 0.0053 0.0051 0.0043 0.0031 0.0016 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 t=1.4时刻的绝对误差0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0002 0.0003 0.0005 0.0006 0.0006 0.0006 0.0005 0.0003 0.0002 0.0000 0.0000 0.0003 0.0007 0.0009 0.0011 0.0011 0.0011 0.0009