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

学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013专业班信计02班学生姓名学号开课时间2015至2016学年第2学期总成绩教师签名2数学与统计学院制开课学院、实验室：数统学院实验时间：2016年6月20日实验项目名称二维波动方程的有限差分法实验项目类型验证演示综合设计其他指导教师曾芳成绩是一．实验目的通过该实验，要求学生掌握求解二维波动方程的有限差分法，并能通过计算机语言编程实现。二．实验内容考虑如下的初值问题：22222222,,0,1,0,1.4,,0sinsin,,,00,,0,1,,0,,,0,1.4uuuxyttxyuxyxyuxyxytuxytxyt(1)1．在第三部分写出问题（1）三层显格式。2．根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。3．取0.1,0.1hh，分别将0.5,1.0,1.4t时刻的数值解画图显示。4.该问题的解析解为,,tcos2sinsinuxytxy，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。三．实验原理、方法（算法）、步骤网格划分0.1,0.1hh，故11.410,140NMh，,,,0,1,,10ijxihyjhij，ktk，0,1,,140k。在内网点,,ijkxyt，利用二阶中心差商，对（1）建立差分格式：11,,,1,,1,,1,,1222222kkkkkkkkkijijijijijijijijijuuuuuuuuuhh（2）整理得到：1221,1,1,,1,1,,24kkkkkkkijijijijijijijuruuuuruu（3）3其中，,1,2,,9,1,2,,139ijk，网比0.1rh，局部截断误差为22oh。考虑边界条件,,0,,,0,1.4uxytxyt，差分格式为：0,00,N,0,0,0,1,,140kkkkNNNuuuuk（4）考虑初始条件,,0sinsinuxyxy，差分格式为：0,sinsinsinsin,,0,1,,10ijijuxyihjhij（5）考虑初始条件2,,00,,0,1tuxyxy，利用二阶差商近似：11,,0,,0,1,,102ijijuuij（6）设0k时刻的点为内点，则满足差分格式（2），代入上式得到：120000201,1,1,,1,1,,24ijijijijijijijuruuuuruu（7）将（6）得到的结果11,,ijijuu代入（7）中，整理得到：12000020,1,1,,1,1,1122ijijijijijijuruuuuru（8）综上（2）、（4）、（5）、（8）得到三层显格式的差分格式为：1221,1,1,,1,1,,0,00,N,0,0,12000,1,1,,124,1,2,,9,1,2,,1390,0,1,,140sinsinsinsin,,0,1,,1012kkkkkkkijijijijijijijkkkkNNNijijijijijijiuruuuuruuijkuuuukuxyihjhijuruuuu020,1,12,,0,1,,10jijruij（9）其中0.1rh，局部截断误差为22oh。四．实验环境（所用软件、硬件等）及实验数据文件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);%空间网格剖分4uu=cos(sqrt(2)*pi*t).*sin(pi*x).*sin(pi*y);%精确解计算%第一层网点计算u=sin(pi*x).*sin(pi*y);%初始条件u1=u(:,:,1);%因为此时得到的u为11x11x141，故只取第一层%第二层网点计算fori=2:10forj=2:10u(i,j,2)=0.5*r^2*(u(i+1,j,1)+u(i-1,j,1)+u(i,j+1,1)+u(i,j-1,1))+(1-2*r^2)*u(i,j,1);u(11,:,2)=0;u(:,11,2)=0;endendu2=u(:,:,2);%第3-141层网点计算fork=2:140fori=2:10forj=2:10u(i,j,k+1)=r^2*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k))+(2-4*r^2)*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时刻的绝对误差矩阵11x115x0=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时刻的精确解');6xlabel('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值');7五．实验结果及实例分析1、0.10.51.01.4t、、、时刻的数值解与精确解图图1t=0.1、0.5时刻的数值解、精确解图2t=1.0、1.4时刻的数值解、精确解注：上两图为四个时刻的数值解与精确解，10.12rpphp代表维数，本文，三层显格式达二阶收敛，不难看出，收敛效果很好，符合理论。下图是四个时刻的绝对误差图像，从图中看出，绝对误差较小，且经过计算得到，收敛阶近似于2，正好符合理论值。82、0.10.51.01.4t、、、时刻的绝对误差图图3四个时刻的绝对误差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.00050.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.00000.00000.00010.00030.00040.00040.00050.00040.00040.00030.00010.00000