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

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

2、00, x,y ,t 0,1.4在第三部分写出问题（1）三层显格式。根据你写出的差分格式，编写有限差分法程序。取 h 0.1,该问题的解析解为u x,y,t示，对数值结果进行简单的讨论。4.0,1.40, x, y20,1(1)将所写程序放到第四部分。0.1h，分别将t 0.5,1.0,1.4时刻的数值解画图显示。cos ；2 tsin xs in y，将四个时刻的数值解的误差画图显网格划分h0.1,0.1h，故 N1 h10, M 口140，人ih, yjtk k，k 0,1,|,140。在内网点x, yj,tk ，利用二阶中心差商，对（k 1_ kk 1k_ kkkkkUi,j2u- 厶i

3、 ,jui,jui 1,j2ui,jUi 1,jui,j12ui,jui,j12h2h2整理得到：k 12kkkkc,2kk1ui,jrui 1,j ui 1,jui,j 1ui,j 124rui,ju,j1）建立差分格式:jh, i,j 0,1,|,10,三.实验原理、方法（算法）、步骤其中，i,j 1,2,川,9, k 1,2,卅,139，网比考虑边界条件u x,y,t 0, x,y,th01，局部截断误差为00,1.4，差分格式为：h2。考虑初始条件u0 ui,jsin考虑初始条件kkkU0,0U0,NUN,0x,y,0X sinx, y,0kUN ,Nsin xsinyj sin0,

4、x,yih0,k 0,1|,140差分格式为:sin jh ,i,j0,1,卅,1020,1 ，利用二阶差商近似:(4)(5)1ui,j21Ui,j0,i, j0,1，川,10(6)设k 0时刻的点为内点，贝U满足差分格式（2），代入上式得到:ui 1i,j1 2 0 0 0 0 2 0ui,j rui 1,j ui 1,j u,j 1 ui,j 12 4r q,j(7)将（6）得到的结果u1,ju,；代入（7）中，整理得到:综上（2）、0Ui,j1Ui,j其中rh1Ui,j(4)、k 1 ui,jsin1 2-r20.11 2-r2(5)、(8)得到三层显格式的差分格式为:2 kkrui 1

5、,j uii, jkku0,0u0,Nxi sin yj0 0 0 0 2 ui 1,j ui 1,j ui,j 1 ui,j 11 2r0 0Ui 1,j Ui 1,j0Ui,j(8)kk2 kk 1ui,j 1ui,j 12 4r ui,j ui,j,9,kkUN ,Nih1,j1,2，川,9,k 1,2,川,1390,k 0,1,|,140sin jh ,i,j 0,1,川,1011 2r2 u0j,i, j 0,1,|,10kUN,0sin(9)0Ui,j0Ui,j局部截断误差为h2。四实验环境（所用软件、硬件等）及实验数据文件Matlab%二维波动方程数值计算（关键：怎么运用i,j,

6、k 三个指标建立循环）clc;如以将代码换成函数m文件h=;tau=*h;%定义步长r=tau/h;% 网比:h:1,0:h:1,0:tau:;%空间网格剖分uu=cos(sqrt(2)*pi*t).*si n(pi*x).*si n(pi*y);%精确解计算%第一层网点计算u=si n( pi*x).*si n(pi*y);%初始条件u仁u(:,:,1);% 因为此时得到的u为11x11x141，故只取第一层%第二层网点计算for i=2:10for j=2:10u(i,j,2)=*rA2*(u(i+1,j,1)+u(i-1,j,1)+u(i,j+1,1)+u(i,j-1,1)+(1-2W2

7、)*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: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%果分%与作%B%wucha=abs(u-uu);% 求绝对误差矩阵 11x11x141wucha仁wucha(:,:,11);% 计算t=时刻的绝对

8、误差矩阵11x11 wucha2=wucha(:,:,51);% 计算t=时刻的绝对误差矩阵11x11 wucha3=wucha(:,:,101);% 计算t=时刻的绝对误差矩阵 11x11 wucha4=wucha(:,:,141);% 计算t=时刻的绝对误差矩阵 11x11 x0=0:h:1;y0=0:h:1;%误差%乍t=时刻的绝对误差图subplot(2,2,1);mesh(xO,yO,wucha1);title('t=时刻的绝对误差');xlabel('x 变量');ylabel('y 变量');zlabel(' %乍t=时刻的

9、绝对误差图subplot(2,2,2);mesh(x0,y0,wucha2);title('t=时刻的绝对误差');xlabel('x 变量');ylabel('y 变量');zlabel('%乍t=时刻的绝对误差图subplot(2,2,3);mesh(x0,y0,wucha3);title('t=时刻的绝对误差');xlabel('x 变量');ylabel('y 变量');zlabel(' %乍t=时刻的绝对误差图subplot(2,2,4);mesh(x0,y0,wucha4

10、);title('t=时刻的绝对误差');xlabel('x 变量');ylabel('y 变量');zlabel('绝对误差值');绝对误差值');绝对误差值');绝对误差值');%四个时刻数值解%精确解 %乍t=、时刻的数值解与精确解subplot(2,2,1);mesh(x0,y0,u(:,:,11);% 作 t=时刻的数值解 title('t=时刻的数值解');xlabel('x 变量');ylabel('y 变量');zlabel('u值&#

11、39;);<)%subplot(2,2,2);mesh(x0,y0,uu(:,:,11);% 作 t=时刻的精确解 title('t= 时刻的精确解');'u%作 t=时刻的数值解与精确解 subplot(2,2,3);mesh(x0,y0,u(:,:,51);% 作 t=时刻的数值解 title('t=时刻的数值解');xlabel('x 变量');ylabel('y变量');zlabel('u值'）；subplot(2,2,4);mesh(x0,y0,uu(:,:,51);% 作 t=时刻的精确解

12、 title('t= 时刻的精确解');xlabel('x 变量');ylabel('y变量');zlabel('u值'）；%制%f%!?%)%作t=、时刻的数值解与精确解subplot(2,2,1);mesh(x0,y0,u(:,:,101);%作t=时刻的数值解title('t=时刻的数值解');xlabel('x 变量');ylabel('y subplot(2,2,2);mesh(x0,y0,uu(:,:,101);% title('t=时刻的精确解');xlabel

13、('x 变量');ylabel('y变量'）;zlabel（'u值'）作t=时刻的精确解变量'）;zlabel（'u值'）%作t=时刻的数值解与精确解subplot(2,2,3); mesh(x0,y0,u(:,:,141);% title('t=时刻的数值解');xlabel('x 变量');ylabel('y subplot(2,2,4);mesh(x0,y0,uu(:,:,141);% title('t=时刻的精确解');xlabel('x 变量');ylabel('y作t=时刻的数值解变量'）;zlabel（'u值'）作t=时刻的精确解变量'）;zlabel（'u值'）五实验结果及实例分析1、t 0.10.51.01.4时刻的数值解与精确解图图1 t=、时刻的数值解、精确解图2 t=、时刻的数值解、精确解注：上两图为四个时刻的数值解与精确解,0.11p代表维数，本文p 2，三层显格式达二阶收敛，不难看出，收敛效果很好，符合理论。下图是四个时刻的