最新一维波动方程的有限差分法.docx
《最新一维波动方程的有限差分法.docx》由会员分享，可在线阅读，更多相关《最新一维波动方程的有限差分法.docx（11页珍藏版）》请在冰点文库上搜索。
最新一维波动方程的有限差分法
一维波动方程的有限差分法
学生实 验 报告
实验课程名称 偏微分方程数值解
开课实验室 数统学院
学院 数统年级 2013专业班信计02班
学生姓 名 学号
开课时 间 2０15至 2016学年第 2 学期
总成 绩
教师签名
数学与统计学院制
开课学院、实验室:
数统学院 实验时间:
２016年６月２0日
实验项目
名 称
一维波动方程的有限差分法
实验项目类型
验证
演示
综合
设计
其他
指导教师
曾芳
成绩
是
一．实验目的
通过该实验，要求学生掌握求解一维波动方程的有限差分法,并能通过计算机语言编程实现。
二.实验内容
考虑如下的初值问题：
（１）
1．在第三部分写出问题
（1）三层显格式。
２.根据你写出的差分格式,编写有限差分法程序．将所写程序放到第四部分。
3.取
，分别将
时刻的数值解画图显示。
4。
该问题的解析解为
将四个时刻的数值解的误差画图显示,对数值结果进行简单的讨论。
三.实验原理、方法（算法）、步骤
1、三层显格式建立
由于题中
取
，故令网比
，
在
内网个点处，利用二阶中心差商得到如下格式:
（２）
略去误差项得到：
（３）
其中
，局部截断误差为
。
对于初始条件
建立差分格式为:
（4）
对于初始条件
利用中心差商,建立差分格式为:
（5）
对于边界条件
建立差分格式为:
（６）
将差分格式延拓使
为内点，代入（３）得到的式子再与（5）联立消去
后整理得到:
（７）
综上（3）、（4）、（6）、（7）得到三层显格式如下:
（局部截断误差为
）
（8）
其中
。
四．实验环境（所用软件、硬件等）及实验数据文件
Ｍaｔlab
三层显格式程序如下：
％一维波动方程，三层显格式求解法
h＝0.１；tau=０.1＊h;
ｒ=tau/h；N=1/h;Ｍ=２/tａu;
x＝0:
h:
１;t=0：
tａｕ:
2；
u=sin（pｉ*x）;％计算t=0时刻的u值
u（1,１1）=0;
forｊ=2：
N
u（2，j）＝0。
5＊r＾2＊u（1,j+1）+（１—r^2）*ｕ（1,j）+０。
5*r^2*u（１，j－1）;
end
%定义x=０边界上的数值
fork=1：
M+1
u（k,1）＝0；
eｎｄ
%定义x=１边界上的数值
fork=1：
M+1
u（ｋ,Ｎ+1）＝0;
end
%迭代计算开始,差分格式
for k=２:
Ｍ
for ｊ=2：
N
u（k+1,j）=r^2*u（ｋ,ｊ+1）+2*（1—r^2）*u（ｋ,j）＋r＾2＊ｕ（ｋ，j—１）—u（ｋ—1,ｊ）;
ｅnd
enｄ
u（201,:
）=zeroｓ（１，11）;
％计算k=2０1行的数值解
u2（２0１,１１）=0；
fｏrj=２:
N
ｕ２（201,ｊ）=ｒ^2＊u（200,ｊ＋1）+２*（1-ｒ^2）＊u（200,j）＋ｒ^2＊u（200,j—1）—u（199，ｊ）；
eｎｄ
ｕ=u+u2;
u=ｒot９0（ｕ,２）;％将矩阵u旋转180度赋值于u
%作出图像
[ｘ，t］=meshgrid（0:
０.1:
1,0:
０.01:
2）;％划分网格
%作出数值解的函数图像
sｕbplot（2，2，１）;
mesh（x，t,u）；
ｔitlｅ（'u（ｘ,t）数值解的函数图像’）；
xlaｂeｌ（'ｘ变量’）;
ylabel（＇t变量’）；
zlaｂel（’u值’）；
%作出精确解的函数图像
suｂplot（２，2，2）；
u１=cos（pi*t）。
*sin（pi*x）；
meｓh（x，t,ｕ１）;
tｉtｌｅ（’u（ｘ,ｔ）精确解的函数图像'）；
ｘlaｂeｌ（’ｘ变量’）；
yｌabel（'t变量＇）；
zlaｂｅl（'u值＇）;
％作出t=0.５,1．０,1．5,2。
０时刻的绝对误差图像
subpｌot（２,２，3）；
wucha=aｂs（u—u1）；
x＝0:
h:
1;
ｐｌoｔ（ｘ，wuｃhａ（5１，:
），'ｇ＊-’）；
hｏld oｎ
griｄｏn
pｌot（ｘ,wuｃｈa（101,:
）,’ro-'）;
hｏｌdon
ｐlot（x,wuｃhａ（151,:
）,'ｋs—'）;
ｈｏldon
plｏt（x,ｗｕｃha（20１,:
）,’mp-＇）;
tｉtle（'t=0。
５,1。
0,1.5,2.0时刻的绝对误差函数图像'）;
xlａbｅl（’x变量'）;ylａbel（'绝对误差值'）;lｅgｅｎｄ（’ｔ=０.5’,'t=1。
０’,'t＝１.5＇,'t＝２。
0’）；
%作出ｔ=0。
５，1．0，1。
5, 2。
０时刻的数值解函数图像
suｂploｔ（２,２,4）;
x=０:
h:
1;
pｌot（x,u（51，:
），’g*—'）;
hｏldon
griｄon
ploｔ（ｘ,u（１01，：
），’ro—’）;
hoｌdｏｎ
ｐlot（x,ｕ（15１,:
），'ｋs-'）；
hｏldｏｎ
ｐloｔ（ｘ,ｕ（201,：
），’mｐ-’）;
ｔitｌe（’ｔ＝0.5,1.0,１.5， 2。
0时刻的数值解函数图像’）；
xlａbｅｌ（'x变量’）；ylabel（'u值＇）;lｅgenｄ（’t＝0。
5’,'t=1。
0'，’t＝1.5'，＇t=2。
0’）;
%当然也可以作出u（x，ｔ）绝对误差的函数图像
%mesh（x,t，wucha）；
％tｉｔｌｅ（'u（ｘ,t）绝对误差的函数图像'）;
％xlａbｅl（’x变量'）;
%ylaｂel（'ｔ变量'）;
%ｚlabｅl（’绝对误差值＇）;
五．实验结果及实例分析
１、u（x,t）在ｔ＝0。
5,1.0,1.5,2．０时刻的数值解、精确解以及绝对误差
表1u（x,ｔ）在t=0。
5，1．0,1。
５,2.０时刻的数值解
时刻ｔ
t=0。
5，1。
0,1。
5,2.0时刻的数值解
ｔ=0。
5
0
—0。
0059
—0。
0１13
—0。
0１55
—０.0182
—０.0１9２
—0.0182
-0。
0155
—0.01１３
—0。
0059
0
t=１。
0
０
—0.3０90
—0。
５8７7
-０。
8090
—０。
9510
－0。
99９9
—0。
９510
-0。
80９0
－０.58７7
-0。
3090
0
t=１.5
0
0。
0020
0。
00３8
0．０05２
0.0061
0.０064
0.０061
0.0０52
0.０0３8
0。
00２0
0
t=2。
0
0
0。
30９０
0。
5878
0．8０90
0．９５11
1。
０００0
0.9511
0.8090
0.5878
0。
3090
0
表２u（x,t）在ｔ=0。
５，1。
0,1。
5，2。
0时刻的精确解
时刻ｔ
t=０．5，1。
0,１。
５，2.０时刻的精确解
t=0.5
0
0。
0000
0。
0０00
0。
0００0
０.０00０
0。
0０00
０。
００0０
0。
000０
0.0０0０
0。
０000
0
t=1。
０
0
－0。
30９０
—０。
5８7８
-0．８０9０
—0.9511
-1。
0000
—０。
９511
—０．809０
—0.５878
－0。
3０90
０
t=1.５
0
0。
０000
0。
０00０
０.0０00
0。
0000
０。
0０００
0.0000
０.０000
0.00０0
０.0000
０
t=２.0
0
0.3090
0。
58７8
０。
80９0
0．951１
1。
０0０0
0.95１1
０.8090
０。
58７8
0.3090
０
表3ｕ（x,t）在ｔ=０.5，1.0，1。
5，2．0时刻的绝对误差
时刻t
t=０。
5,１.０,1。
5，２。
0时刻的绝对误差
t=0.５
0
０。
0059
０。
01１３
0．01５５
0.０１82
0．0１9２
0.0182
0.０155
0。
0113
0。
0０５9
0
t=1.0
0
0.0000
0。
00０0
0．0001
0。
00０1
０。
0001
0.0001
0.000１
0.0000
0.000０
0
t=1。
5
0
0．00２０
０。
０0３８
0。
００52
0．0061
0.0064
0。
00６1
０。
0052
0.0０38
０.00２0
０
t=2。
0
０
0。
0０００
0．0000
0．０００0
0.0０00
0．0０00
0.0000
0。
0000
０。
000０
0。
0０0０
0
说明:
在t=0。
5时刻的绝对误差最大，t=1。
５时刻次之,t=１与ｔ=2时刻的绝对误差均较小,由于
，该格式稳定,由数值计算得到的矩阵不难看出,数值解符合理论解.
2、ｕ（x,t）在t=0.5,1。
0,1。
5，2。
0时刻的数值解、绝对误差函数图像
图1数值解、精确解以及绝对误差函数图像
说明：
上两图为函数的数值解与精确解,下两图为ｔ=0。
5,１.0,1。
5,2．0时刻的数值解、绝对误差函数图像，符合理论解。
教师签名
年 月 日
升级会员 