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

学生实验报告

实验课程名称偏微分方程数值解

开课实验室数统学院

学院数统年级2013专业班信计02班

学生姓名学号

开课时间2015至2016学年第2学期

总成绩教师签名

数学与统计学院制

开课学院、实验室：数统学院实验时间：2016年6月20日
实验项目
名称一维波动方程的有限差分法实验项目类型验证演示综合设计其他指导教师曾芳成绩是一．实验目的
通过该实验，要求学生掌握求解一维波动方程的有限差分法，并能通过计算机语言编程实现。

二．实验内容
考虑如下的初值问题：

(1)
1．在第三部分写出问题（1）三层显格式。
2．根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。
3．取，分别将时刻的数值解画图显示。
4.该问题的解析解为，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。

三．实验原理、方法（算法）、步骤

1、三层显格式建立
由于题中，，取，故令网比，，，在内网个点处，利用二阶中心差商得到如下格式：
（2）
略去误差项得到：
（3）
其中，局部截断误差为。
对于初始条件，建立差分格式为：
（4）
对于初始条件，利用中心差商，建立差分格式为：
（5）
对于边界条件，建立差分格式为：
（6）
将差分格式延拓使为内点，代入（3）得到的式子再与（5）联立消去后整理得到：
（7）
综上（3）、（4）、（6）、（7）得到三层显格式如下：（局部截断误差为）
（8）
其中。

四．实验环境（所用软件、硬件等）及实验数据文件
Matlab
三层显格式程序如下：
%一维波动方程，三层显格式求解法
h=0.1;tau=0.1*h;
r=tau/h;N=1/h;M=2/tau;
x=0:h:1;t=0:tau:2;
u=sin(pi*x);%计算t=0时刻的u值
u(1,11)=0;
forj=2:N
u(2,j)=0.5*r^2*u(1,j+1)+(1-r^2)*u(1,j)+0.5*r^2*u(1,j-1);
end
%定义x=0边界上的数值
fork=1:M+1
u(k,1)=0;
end
%定义x=1边界上的数值
fork=1:M+1
u(k,N+1)=0;
end
%迭代计算开始,差分格式
fork=2:M
forj=2:N
u(k+1,j)=r^2*u(k,j+1)+2*(1-r^2)*u(k,j)+r^2*u(k,j-1)-u(k-1,j);
end
end
u(201,:)=zeros(1,11);
%计算k=201行的数值解
u2(201,11)=0;
forj=2:N
u2(201,j)=r^2*u(200,j+1)+2*(1-r^2)*u(200,j)+r^2*u(200,j-1)-u(199,j);
end
u=u+u2;
u=rot90(u,2);%将矩阵u旋转180度赋值于u
%作出图像
[x,t]=meshgrid(0:0.1:1,0:0.01:2);%划分网格
%作出数值解的函数图像
subplot(2,2,1);
mesh(x,t,u);
title('u(x,t)数值解的函数图像');
xlabel('x变量');
ylabel('t变量');
zlabel('u值');
%作出精确解的函数图像
subplot(2,2,2);
u1=cos(pi*t).*sin(pi*x);
mesh(x,t,u1);
title('u(x,t)精确解的函数图像');
xlabel('x变量');
ylabel('t变量');
zlabel('u值');
%作出t=0.5，1.0，1.5,2.0时刻的绝对误差图像
subplot(2,2,3);
wucha=abs(u-u1);
x=0:h:1;
plot(x,wucha(51,:),'g*-');
holdon
gridon
plot(x,wucha(101,:),'ro-');
holdon
plot(x,wucha(151,:),'ks-');
holdon
plot(x,wucha(201,:),'mp-');
title('t=0.5，1.0，1.5,2.0时刻的绝对误差函数图像');
xlabel('x变量');ylabel('绝对误差值');legend('t=0.5','t=1.0','t=1.5','t=2.0');
%作出t=0.5，1.0，1.5,2.0时刻的数值解函数图像
subplot(2,2,4);
x=0:h:1;
plot(x,u(51,:),'g*-');
holdon
gridon
plot(x,u(101,:),'ro-');
holdon
plot(x,u(151,:),'ks-');
holdon
plot(x,u(201,:),'m