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

1、学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013专业班信计02班学生姓名学号开课时间2015至2016学年第2学期总 成 绩教师签名.数学与统计学院制开课学院 、实验室： 数统学院实验时间 ：2016 年 6月20日实验项目实验项目类型二维波动方程的有限差分法名称验证演示综合设计其他指导教师曾芳成绩是一实验目的通过该实验 ，要求学生掌握求解二维波动方程的有限差分法 ，并能通过计算机语言编程实现。二实验内容考虑如下的初值问题 ：2u2u2u2, x, y0,1.40,1 , tt 2x2y2ux, y,0sinx siny,u x, y,00, x, y0,12(1

2、)tux, y,t0, x, y, t0,1.41在第三部分写出问题 （1）三层显格式 。2根据你写出的差分格式 ，编写有限差分法程序 。 将所写程序放到第四部分。3取 h0.1,0.1h ，分别将 t0.5,1.0,1.4 时刻的数值解画图显示。4. 该问题的解析解为 u x, y,t cos 2 t sin x sin y ，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论 。三实验原理 、方法（算法）、步骤.专业资料 .网格划分h0.1,0.1h，故110, M1.4，N140hxi ih , y jjh , i , j0,1,10 ， tkk， k0,1,140 。在内网点x

3、i , yj , tk，利用二阶中心差商，对（1）建立差分格式 ：uik, j 12uik, juik, j 1uik 1, j2uik, juik 1, juik, j 1 2uik, juik, j 1（2）2h2h2整理得到 ：uik, j 1r 2 uik 1, juik 1, juik, j 1 uik, j 12 4r 2 uik, juik, j 1（3）其中， i, j 1,2, ,9, k1,2,139 ，网比 rh0.1，局部截断误差为 o2h2 。考虑边界条件 u x, y,t0, x, y, t0,1.4 ，差分格式为 ：u0,0ku0,NkuNk,0 uNk,N0,k

4、0,1, ,140考虑初始条件 ux, y,0sinx siny ，差分格式为 ：u0jsinxisinyjsinih sinjh ,i , j 0,1, ,10i,考虑初始条件 utx, y,00,x, y0,12 ，利用二阶差商近似 ：ui1, jui ,1j0, i , j0,1,102设 k0 时刻的点为内点 ，则满足差分格式 （2），代入上式得到 ：ui1, jr 2 ui0 1, jui0 1, jui0, j 1 ui0, j 12 4r 2 ui0, j ui ,1j将（ 6）得到的结果 ui1, jui ,1j 代入（7）中，整理得到 ：ui1, j1 r 2ui01,j u

5、i01, j ui0, j 1 ui0, j 11 2r 2 ui0, j2综上（2）、（ 4）、（ 5）、（ 8）得到三层显格式的差分格式为：uik, j 1r 2 uik 1, juik1, juik, j 1uik, j 124r 2 uik, juik, j 1i , j1,2,9, k1,2,139u0,0ku0,NkuNk ,0uNk ,N0, k0,1,140ui0, jsinxisiny jsinihsinjh ,i , j0,1,10ui1, j1 r 2ui01, jui01, jui0, j1ui0, j112r2 ui0, j ,i , j0,1, ,102（4）（5）

6、（6）（7）（8）（9）.专业资料 .其中 r0.1 ，局部截断误差为 o2h2 。h四实验环境 （所用软件 、硬件等）及实验数据文件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

7、(:,:,1);% 因为此时得到的u 为 11x11x141 ，故只取第一层%第二层网点计算for i=2:10for j=2:10u(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:140for i=2:10for 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

8、-4*r2)*u(i,j,k)-u(i,j,k-1);u(11,:,k+1)=0;u(:,11,k+1)=0;endendend%结果分析与作图%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

9、: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);

10、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时刻的数

11、值解 );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,

12、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);%

13、 作 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 变量 );ylabe

14、l(y 变量 );zlabel(u 值);.专业资料 .五实验结果及实例分析1、 t0.1、0.5、1.0、1.4 时刻的数值解与精确解图图 1 t=0.1 、0.5 时刻的数值解 、精确解.专业资料 .图 2 t=1.0 、1.4 时刻的数值解 、精确解注：上两图为四个时刻的数值解与精确解， r0.11p代表维数 ，本文 p 2hp，三层显格式达二阶收敛 ，不难看出 ，收敛效果很好 ，符合理论 。 下图是四个时刻的绝对误差图像 ，从图中看出 ，绝对误差较小 ，且经过计算得到 ，收敛阶近似于 2，正好符合理论值。2、 t0.1、0.5、1.0、1.4 时刻的绝对误差图.专业资料 .图 3 四个

15、时刻的绝对误差3、四个时刻 （ t=0.1 、0.5 、1.0 、 1.4 ）的绝对误差表t=0.1时刻的绝对误差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00010.00010.00020.00020.00020.00020.00020.00010.00010.00000.00000.00010.00030.00040.00040.00050.00040.00040.00030.00010.00000.00000.00020.00040.00050.00060.00060.00060.0

16、0050.00040.00020.00000.00000.00020.00040.00060.00070.00070.00070.00060.00040.00020.00000.00000.00020.00050.00060.00070.00080.00070.00060.00050.00020.00000.00000.00020.00040.00060.00070.00070.00070.00060.00040.00020.00000.00000.00020.00040.00050.00060.00060.00060.00050.00040.00020.0000.专业资料 .0.00000.

17、00010.00030.00040.00040.00050.00040.00040.00030.00010.00000.00000.00010.00010.00020.00020.00020.00020.00020.00010.00010.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000t=0.5时刻的绝对误差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00070.00130.00180.00210.

18、00220.00210.00180.00130.00070.00000.00000.00130.00250.00340.00400.00420.00400.00340.00250.00130.00000.00000.00180.00340.00470.00550.00580.00550.00470.00340.00180.00000.00000.00210.00400.00550.00650.00680.00650.00550.00400.00210.00000.00000.00220.00420.00580.00680.00710.00680.00580.00420.00220.00000.

19、00000.00210.00400.00550.00650.00680.00650.00550.00400.00210.00000.00000.00180.00340.00470.00550.00580.00550.00470.00340.00180.00000.00000.00130.00250.00340.00400.00420.00400.00340.00250.00130.00000.00000.00070.00130.00180.00210.00220.00210.00180.00130.00070.00000.00000.00000.00000.00000.00000.00000.

20、00000.00000.00000.00000.0000t=1.0时刻的绝对误差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00160.00310.00430.00510.00530.00510.00430.00310.00160.00000.00000.00310.00590.00820.00960.01010.00960.00820.00590.00310.00000.00000.00430.00820.01130.01320.01390.01320.01130.00820.00430.

21、00000.00000.00510.00960.01320.01560.01640.01560.01320.00960.00510.00000.00000.00530.01010.01390.01640.01720.01640.01390.01010.00530.00000.00000.00510.00960.01320.01560.01640.01560.01320.00960.00510.00000.00000.00430.00820.01130.01320.01390.01320.01130.00820.00430.00000.00000.00310.00590.00820.00960.01010.00960.00820.00590.00310.00000.00000.00160.00310.00430.00510.00530.00510.00430.00310.00160.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000t=1.4时刻的绝对误差0.00000.00000.00000.00000.0000