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

1、学生实验报告实验课程名称偏微分方程数值解开课实验室 数统学院学院 数统年级 20院专业班 信计02班学生姓名 学 号开课时 间 2015至2016学年第 2 学期总成绩教师签名数学与统计学院制x,y,t=meshgrid(0:h:1,0:h:1,0:tau:1.4);%空间网格剖分开课学院、实验室： 数统学院实验时间：2016年6月20日实验项目名 称二维波动方程的有限差分法实验项目类型验证演示综合设计其他指导教师目力成 绩是一.实验目的通过该实验，要求学生掌握求解二维波动方程的有限差分法， 并能通过计算机语言编程 实现。.实验内容考虑如下的初值问题:-2_2_2鬟=劈十粤(x,y)wc =(

2、0,1)2,t10,1.4):t二 x: y22«u (x, y,0 产sinxsin y y,一u (x, y,0 )= 0,(x, y 产(0,1)x,y,t =0, x,y r,t10,1.41在第三部分写出问题（1）三层显格式。2.3.4.根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。取h =0.1,7 =0.1h ,分别将t =0.5,1.0,1.4时刻的数值解画图显示。该问题的解析解为u（x,y,t ） = cosJ2ntsinnxsinny ,将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。.实验原理、方法（算法）、步骤网格划分h =0.

3、1,7 =0.1h，故N =1.4htk =kT , k =0,1,111,140。在内网点(K, yj,tk )= 10, M = 140, Xi =ih,yj = jh, i, j =0,1川，10 ,利用二阶中心差商，对（1）建立差分格式:k 1 八 kUi,j-2ui,j2k -1kk kkk k' ui,j =ui1,j-2ui,j ' ui 4,j . ui,j 1 - 2ui,j ' ui,j 4h2h2(2)整理得到:k 12 kkui,j =r ui u ' IUM Uikj42-4r2 u：j y7(3)其中，i, j =12|,9,k=12

4、|,139 ,网比 r=?=0.1,局部截断误差为 o（d + h2）。 h考虑边界条件u x,y,t =0, x,y）三1"t . 0,1.41,差分格式为：_ _k _ _ k_ _ k_ _ ku0,0 = u0,N = uN,0 = uN ,N = 0,k =0,1,H|,140(4)考虑初始条件u（x,y,0）= sinnxsinny ,差分格式为:0ui,j =sin 二 xi sin r： yj =sin 二 ih sin 二 jh ,i, j = 0,1,l|l ,10(5)2考虑初始条件ut（x,y,0 ）=0,（x,y产（0,1）,利用二阶差冏近似:1 jUi j

5、 - Uj j,j,j =0,i, j =0,1,111,102.设k=0时刻的点为内点，则满足差分格式(2),代入上式得到:u1,j =r2(u*j +u0,j +Wj+u0j)+(24r2 u0j uj(6)将（6）得到的结果u：j =5：代入（7）中，整理得到:1ui,j综上（2）、(4)、k 1uui,j12000020=2rUi 1,jujui,j 1 Ui,j1-2rui,j(5)、(8)得到三层显格式的差分格式为：2kkkk2k krui 1,j uj "j 1 U,j2-4r uJ-u5(8)0 u ui,ju1,ji,j =1,2,|l|,9,k =1,2,小，13

6、9k k k ku0,0 = u0,N =uN,0 =uN,N - 0, k - 0,1,| I I ,140二sin - xi sin 二 yj l =sin 二ih sin 二 jh ,i, j =0,1,|l|,10= 2r2(u0+j +u0,j +u0j + +ui0j)+(12r2 )u0j,i, j =0,1川,10(9)其中=三=0.1,局部截断误差为o（T2+h2）。 h四.实验环境（所用软件、硬件等）及实验数据文件Matlab%二维波动方程数值计算（关键：怎么运用i,j,k 三个指标建立循环）clc;%可以将代码换成函数m文件h=0.1;tau=0.1*h;% 定义步长 r

7、=tau/h;% 网比uu=cos(sqrt(2)*pi*t).*sin(pi*x).*sin(pi*y);%精确解计算%第一层网点计算u=sin(pi*x).*sin(pi*y);% 初始条件u1=u(:,:,1);%因为此时得到的u为11x11x141,故只取第一层%第二层网点计算fo门=2:10for j=2:10u(i,j,2)=0.5*rA2*(u(i+1,j,1)+u(i-1,j,1)+u(i,j+1,1)+u(i,j-1,1)+(1-2*rA2)*u(i,j, 1)；u(11,:,2)=0;u(:,11,2)=0;endendu2=u(:,:,2);%S 3-141层网点计算fo

8、r k=2:140for i=2:10for j=2:10u(i,j,k+1)=rA2*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k)+(2-4*rA2)*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

9、);% 计算t=1.0时刻的绝对误差矩阵11x11wucha4=wucha(:,:,141);% 计算 t=1.4 时亥1I的绝对误差矢口阵 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）;ti

10、tle（'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=i.4时刻的绝对误差图 subplot（2,2,4）;mesh（x0,y0,wucha4）

11、;title（'t=1.4时亥1J的绝对误差'）;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（

12、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

13、,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

14、（'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（'

15、;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、t =0.1、051.01.4时刻的数值解与精确解图图1 t=0.1、0.5时刻的数值解、精确解图2 t=1.0、1.4时刻的数值解、精确解注：上两图为四个时刻的数值解与精确解,1r=0.1<j=（

16、p代表维数）,本文p = 2,三层显格式达二阶收敛，不难看出，收敛效果很好，符合理论。下图是四个时刻的绝对误差图像，从图中看出，绝对误差较小，且经过计算得到，收敛阶近似于2,正好符合理论值2、t=0.1、051.R 1.4时刻的绝对误差图t=1。划的拾"一蝴弱的施好总NQ1时用的也甘心及M 5m的箔灯我把n W0.02作小序.0.1?10006 .0 .-3,1口图3四个时刻的绝对误差3、四个时刻（t=0.1、0.5、1.0、1.4）的绝对误差表t=0.1时刻的绝对误差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00

17、000.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.00050.00040.00020.00000.00000.00020.00040.00060.00070.00070.00070.00060.00040.00020.00000.00000.00020.00050.00060.00

18、070.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.00000.00000.00010.00030.00040.00040.00050.00040.00040.00030.00010.00000.00000.00010.00010.00020.00020.00020.00020.00020.00010.00010.00

19、000.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.00220.00210.00180.00130.00070.00000.00000.00130.00250.00340.00400.00420.00400.00340.00250.00130.00000.00000.00180.00340.00

20、470.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.00000.00210.00400.00550.00650.00680.00650.00550.00400.00210.00000.00000.00180.00340.00470.00550.00580.00550.00470.00340.00

21、180.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.00000.00000.00000.00000.0000t=1.0时刻的绝对误差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00160.00

22、310.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.00000.00000.00510.00960.01320.01560.01640.01560.01320.00960.00510.00000.00000.00530.01010.01390.01640.01720.01640.01390.01

23、010.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.00