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

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

2、1(1)在第三部分写出问题（1）三层显格式。根据你写出的差分格式，编写有限差分法程序。取 h 0.1,该问题的解析解为u x,y,t示，对数值结果进行简单的讨论。4.将所写程序放到第四部分。0.1h，分别将t 0.5,1.0,1.4时刻的数值解画图显示。cosJ2 tsin xs in y，将四个时刻的数值解的误差画图显三. 实验原理、方法（算法）、步骤11 4网格划分 h 0.1,0.1h，故 N - 10,M 一 140，Xiih, yjhjh, i,j 0,1,L ,10 ,tk k ，k0,1丄,140。在内网点Xi,yj,tk ,利用二阶中心差商，对（1 ）建立差分格式:k 1Ui,

3、jkk 12Ui,jUi,jkUi 1,jc k2Ui,jkUi 1,jkkkUi,j 12Ui,jUi,j 12h2h2整理得到：k 1ui,j2kr ui 1,jkui 1,jkUi,j 1kui,j 1CX 2kk 12 4rui,jU,j其中，i,j 1,2,L,9, k 1,2, L,139，网比r h0.1，局部截断误差为考虑边界条件u x,y,t0, x,y,t 0,1.4，差分格式为：kkU0,0U0,NUN,0kUN ,N0,k0,1,L ,140考虑初始条件U x,y,0sin xsin y，差分格式为：u0. Ui, jsinxi sinyjsinih sinjh ,i,

4、j 0,1,L ,10考虑初始条件Ut X, y,020, x, y0,1 ，利用二阶差商近似：1ui,j1ui,j0,i, j0,1,L,1020h2。设k 0时刻的点为内点，贝W足差分格式（2），代入上式得到:120000C/201ui,j rui 1,j ui 1,j u,j 1 ui,j 12 4r ui,j ui,j(7)将（6）得到的结果Sj Ui,1代入（7、中，整理得到:1ui,j1 2 0 0 0 0 .20 2r Ui 1,J Ui 1,J Ui,j 1 Ui,j 11 2r Ui,j综上（2）、(4 )、（5）、（8）得到三层显格式的差分格式为:0Ui,jUi,jk 1u

5、i,jsin2r20.1 ,2 kkkkc/2kr Ui 1,j Ui 1,j Ui,j 1Ui,j 12 4r 嘔i,j 1,2,L ,9,k 1,2,L ,139kkkU0,0U0,NUN,0Xi sin yj sinuN,N 0,k0,1,L ,140ih sin jh ,i, j0,1,L0 0 0ui 1,jUi 1,j Ui,j 1局部截断误差为00“Ic 20 Ui,j 112rUi,j,i, J2 h2。k 1Ui,j,100,1L ,10四. 实验环境(所用软件、硬件等)及实验数据文件Matlab%二维波动方程数值计算(关键：怎么运用i,j,k三个指标建立循环)cic;%可以

6、将代码换成函数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)*p i*t).*si n(pi *x).*si n(p );%第一层网点计算u=si n(pi *x).*si n(pi *y);%初始条件u1=u(:.:.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

7、.j-1.1)+(1-2*r2)*u(i.j.1); u(11.:.2)=0;u(:.11.2)=0;end end u2=u(:,:,2);%第3-141层网点计算 for 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(:,

8、:,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时刻的绝对误差图sub plot(2,2,1);mesh(x0,y0,wucha1);title('t=0.1时刻的绝对误差');xlabel('x 变量');ylabel('y 变量&#

9、39;);zlabel('%作t=0.5时刻的绝对误差图绝对误差值'）;sub plot(2,2,2);mesh(x0,y0,wucha2);title('t=0.5时刻的绝对误差');xlabel('x 变量');ylabel('y 变量');zlabel('%作t=1.0时刻的绝对误差图绝对误差值'）;sub plot(2,2,3);mesh(x0,y0,wucha3);title（'t=1.0时刻的绝对误差'）;xlabel（'x 变量'）;ylabel（'y 变量&#

10、39;）;zlabel（'%作t=1.4时刻的绝对误差图绝对误差值'）;sub plot(2,2,4);mesh(x0,y0,wucha4);title（'t=1.4时刻的绝对误差'）;xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'绝对误差值'）;% 四个时刻数值解、精确 解%作t=0.1、0.5时刻的数值解与精确解 sub plot（2,2,1）;mesh（x0,y0,u（:,:,11）;% 作 t=0.1 时刻的数值解 title（'t=0.1时刻的数值解'）;xl

11、abel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'）; sub plot（2,2,2）;mesh（x0,y0,uu（:,:,11）;% 作 t=0.1 时刻的精确解 title（'t=0.1时刻的精确解'）;xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'）;%作t=0.5时刻的数值解与精确解 sub plot（2,2,3）;mesh（x0,y0,u（:,:,51）;% 作 t=0.5 时刻的数值解 title（

12、9;t=0.5 时刻的数值解'）;xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'）;sub plot（2,2,4）;mesh（x0,y0,uu（:,:,51）;% 作 t=0.5 时刻的精确解 title（'t=0.5时刻的精确解'）;xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'）;% 分别复希 y 粘贴运 行%作t=1.0、1.4时刻的数值解与精确解 sub plot（2,2,1）;mesh（x0

13、,y0,u（:,:,101）;% 作 t=1.0 时刻的数值解 title（'t=1.0 时刻的数值解'）;xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'）; sub plot（2,2,2）;mesh（x0,y0,uu（:,:,101）;% 作 t=1.0 时刻的精确解 title（'t=1.0 时刻的精确解）；xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'）;%作t=1.4时刻的数值解与精确解 su

14、b plot（2,2,3）;mesh（x0,y0,u（:,:,141）;% 作 t=1.4 时刻的数值解 title（'t=1.4时刻的数值解'）;xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'）;sub plot（2,2,4）;mesh（x0,y0,uu（:,:,141）;% 作 t=1.4 时刻的精确解 title（'t=1.4时刻的精确解'）;xlabel（'x 变量'）;ylabel（'y 变量'）;zlabel（'u 值'

15、;）;五. 实验结果及实例分析1、t 0.1、0.51.01.4时刻的数值解与精确解图a丄一v*t+Hfl0107=D4图 1 t=0.1asat04Q0H-ArTHOJ(Siff、0.5时刻的数值解、精确解.1iQi b-1 fi."ri 'bi' *區打穿号<12fl 1 站 i t '.ih HK0204.厶鯉4fl Op：*- 卩諾曲利6J41D4一404DBcert詁图2 t=1.0、1.4时刻的数值解、0.17Pp代表维数，本文P 2注：上两图为四个时刻的数值解与精确解，r -h，三层显格式达二阶收敛，不难看出，收敛效果很好，符合理论。下图是

16、四个时刻的绝对误 差图像，从图中看出，绝对误差较小，且经过计算得到，收敛阶近似于2，正好符合理论值。2、t 0.1、0.51.01.4时刻的绝对误差图f1Du_fltil* 出”口.斧时 ilu&i>ieC'DQ-T；«訥.mon鳥- A亠、JS=.40D*O'«Q»0i*¥图3四个时刻的绝对误差3、四个时刻（t=0.1、0.5、1.0、1.4 ）的绝对误差表t=o.i时刻的绝对误差0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0001

17、0.00010.00020.00020.00020.00020.00020.00010.00010.00000.00010.00030.00040.00040.00050.00040.00040.00030.00010.00000.00020.00040.00050.00060.00060.00060.00050.00040.00020.00000.00020.00040.00060.00070.00070.00070.00060.00040.00020.00000.00020.00050.00060.00070.00080.00070.00060.00050.00020.00000.0002

18、0.00040.00060.00070.00070.00070.00060.00040.00020.00000.00020.00040.00050.00060.00060.00060.00050.00040.00020.00000.00010.00030.00040.00040.00050.00040.00040.00030.00010.00000.00010.00010.00020.00020.00020.00020.00020.00010.00010.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000.0000.

19、0000.0000.0000.0000.0000.0000.0000.0000.00000000000000t=0.5时刻的绝对误差0.0000.0000.0000.0000.0000.0000.0000.0000.0000.0000.000000000000000.0000.0000.0010.0010.0020.0020.0020.0010.0010.0000.000073812183700.0000.0010.0020.0030.0040.0040.0040.0030.0020.0010.000035402045300.0000.0010.0030.0040.0050.0050.0050

20、.0040.0030.0010.000084758574800.0000.0020.0040.0050.0060.0060.0060.0050.0040.0020.000010558550100.0000.0020.0040.0050.0060.0070.0060.0050.0040.0020.000022881882200.0000.0020.0040.0050.0060.0060.0060.0050.0040.0020.000010558550100.0000.0010.0030.0040.0050.0050.0050.0040.0030.0010.000084758574800.0000

21、.0010.0020.0030.0040.0040.0040.0030.0020.0010.000035402045300.0000.0000.0010.0010.0020.0020.0020.0010.0010.0000.000073812183700.0000.0000.0000.0000.0000.0000.0000.0000.0000.0000.00000000000000t=1.0时刻的绝对误差0.0000.0000.0000.0000.0000.0000.0000.0000.0000.0000.000000000000000.0000.0010.0030.0040.0050.005

22、0.0050.0040.0030.0010.000061313131600.0000.0030.0050.0080.0090.0100.0090.0080.0050.0030.000019261629100.0000.0040.0080.0110.0130.0130.0130.0110.0080.0040.000032329232300.0000.0050.0090.0130.0150.0160.0150.0130.0090.0050.000016264626100.0000.0050.0100.0130.0160.0170.0160.0130.0100.0050.00003194249130

23、0.0000.0050.0090.0130.0150.0160.0150.0130.0090.0050.000016264626100.0000.0040.0080.0110.0130.0130.0130.0110.0080.0040.000032329232300.0000.0030.0050.0080.0090.0100.0090.0080.0050.0030.000019261629100.0000.0010.0030.0040.0050.0050.0050.0040.0030.0010.000061313131600.0000.0000.0000.0000.0000.0000.0000.0000