一维波动方程的有限差分法.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;
for j=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边界上的数值
for k=1:M+1
u(k,1)=0;
end
%定义x=1边界上的数值
for k=1:M+1
u(k,N+1)=0;
end
%迭代计算开始,差分格式
for k=2:M
for j=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;
for j=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*-');
hold on
grid on
plot(x,wucha(101,:),'ro-');
hold on
plot(x,wucha(151,:),'ks-');
hold on
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*-');
hold on
grid on
plot(x,u(101,:),'ro-');
hold on
plot(x,u(151,:),'ks-');
hold on
plot(x,u(201,:),'mp-');
title('t=0.5，1.0，1.5, 2.0时刻的数值解函数图像');
xlabel('x变量'); ylabel('u值');legend('t=0.5','t=1.0','t=1.5','t=2.0');
%当然也可以作出u(x,t)绝对误差的函数图像
%mesh(x,t,wucha);
%title('u(x,t)绝对误差的函数图像');
%xlabel('x变量');
%ylabel('t变量');
%zlabel('绝对误差值');
五．实验结果及实例分析

1、u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解、精确解以及绝对误差
表1 u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解
时刻t
t=0.5,1.0,1.5,2.0时刻的数值解
t=0.5
0
-0.0059
-0.0113
-0.0155
-0.0182
-0.0192
-0.0182
-0.0155
-0.0113
-0.0059
0
t=1.0
0
-0.3090
-0.5877
-0.8090
-0.9510
-0.9999
-0.9510
-0.8090
-0.5877
-0.3090
0
t=1.5
0
0.0020
0.0038
0.0052
0.0061
0.0064
0.0061
0.0052
0.0038
0.0020
0
t=2.0
0
0.3090
0.5878
0.8090
0.9511
1.0000
0.9511
0.8090
0.5878
0.3090
0

表2 u(x,t)在t=0.5,1.0,1.5,2.0时刻的精确解
时刻t
t=0.5,1.0,1.5,2.0时刻的精确解
t=0.5
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0
t=1.0
0
-0.3090
-0.5878
-0.8090
-0.9511
-1.0000
-0.9511
-0.8090
-0.5878
-0.3090
0
t=1.5
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0
t=2.0
0
0.3090
0.5878
0.8090
0.9511
1.0000
0.9511
0.8090
0.5878
0.3090
0

表3 u(x,t)在t=0.5,1.0,1.5,2.0时刻的绝对误差
时刻t
t=0.5,1.0,1.5,2.0时刻的绝对误差
t=0.5
0
0.0059
0.0113
0.0155
0.0182
0.0192
0.0182
0.0155
0.0113
0.0059
0
t=1.0
0
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0000
0.0000
0
t=1.5
0
0.0020
0.0038
0.0052
0.0061
0.0064
0.0061
0.0052
0.0038
0.0020
0
t=2.0
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0

说明：在t=0.5时刻的绝对误差最大，t=1.5时刻次之，t=1与t=2时刻的绝对误差均较小，由于，该格式稳定，由数值计算得到的矩阵不难看出，数值解符合理论解。
2、u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解、绝对误差函数图像

图1 数值解、精确解以及绝对误差函数图像

说明：上两图为函数的数值解与精确解，下两图为 t=0.5,1.0,1.5,2.0时刻的数值解、绝对误差函数图像，符合理论解。

教师签名

年 月 日

11