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

学 生 实 验 报 告 实验课程名称 偏微分方程 数统学院 学 院 数 统 年级 2013 专业班 信计02班 学 生 姓 名 学 号 开 课 时 间 2015 至 2016学年第 2 学期 总 成 绩 教师签名 数学与统计学院制 开课学院、实验室： 数统学院 2016年 6月20日 实验项目 名 称 二维波动方程的有限差分法 实验项目类型 验证 演示 综合 设计 其他 指导教师 曾芳 成 绩 是 一．实验目的 通过该实验，要求学生掌握法并能通过计算机语言编程实现。如下的初值问题： (1) 1．在第三部分写出问题（1）三层显格式。 2．根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。 3．取，分别将时刻的数值解画图显示。 4. 该问题的解析解为将时刻的数值解的误差画图显示，的讨论，故，，，。在内网点，利用二阶中心差商，对（1）建立差分格式： （2） 整理得到： （3） 其中，，网比，局部截断误差为。 考虑边界条件，差分格式为： （4） 考虑初始条件，差分格式为： （5） 考虑初始条件，利用二阶差商近似： （6） 设时刻的点为内点，则满足差分格式（2），代入上式得到： （7） 将（6）得到的结果代入（7）中，整理得到： （8） 综上（2）、（4）、（5）、（8）得到三层显格式的差分格式为： （9） 其中，局部截断误差为。 四．实验环境（所用软件、硬件等）及实验数据文件 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(:,:,1);%因为此时得到的u为11x11x141，故只取第一层 %第二层网点计算 for i=2:10 for j=2:10 u(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; end end u2=u(:,:,2); %第3-141层网点计算 for k=2:140 for i=2:10 for j=2:10 u(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; end end end %%%%%%%%%%%%%%%%%%%%%%%%%%结果分析与作图%%%%%%%%%%%%%%%%%%%%%%%%%%%% wucha=abs(u-uu);%求绝对误差矩阵11x11x141 wucha1=wucha(:,:,11);%计算t=0.1时刻的绝对误差矩阵11x11 wucha2=wucha(:,:,51);%计算t=0.5时刻的绝对误差矩阵11x11 wucha3=wucha(:,:,101);%计算t=1.0时刻的绝对误差矩阵11x11 wucha4=wucha(:,:,141);%计算t=1.4时刻的绝对误差矩阵11x11 x0=0:h:1;y0=0:h:1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%误差分析%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %作t=0.1时刻的绝对误差图 subplot(2,2,1); mesh(x0,y0,wucha1); title(t=0.1时刻的绝对误差); xlabel(x变