物理学毕业论文_-_波动方程的数值解法研究.doc
- 文档编号:30816330
- 上传时间:2024-01-30
- 格式:DOC
- 页数:31
- 大小:1.54MB
物理学毕业论文_-_波动方程的数值解法研究.doc
《物理学毕业论文_-_波动方程的数值解法研究.doc》由会员分享,可在线阅读,更多相关《物理学毕业论文_-_波动方程的数值解法研究.doc(31页珍藏版)》请在冰豆网上搜索。
学院
本科毕业论文(设计)
题目波动方程的数值解法研究
院系物理与电子工程学院
专业物理学
姓名
学号
学习年限2007年9月至2011年7月
指导教师
申请学位理学学士学位
二零一一年五月二十三日
波动方程的数值解法研究
摘要:
本文利用有限差分法及其显示差分格式表示了波动方程的差分格式以及初始、边界条件的差分格式,进而讨论了不同初始、边界条件下的波动方程的解的情况,并利用MATLAB数学软件编程绘制出了各种解的图像.为了方便对振动过程的进一步研究,还给出了振动过程的动画演示.对于方程初始、边界条件比较复杂因而无法求得解析解的波动方程,有限差分方法的数值解法给出了一种较为普适而精确的计算方法.
关键词:
波动方程;有限差分法;显式差分格式;数值解
NumericalSolutionofWaveEquation
ABSTRACT:
Inthispaper,usingfinitedifferencemethodandshowdifferenceschemeshowedthatthedifferenceschemeforwaveequationandtheinitialandboundaryconditions,thendiscussedthedifferentinitialandboundaryconditionsofthewaveequationunderthecircumstances,andusingmathematicalsoftwareprogramgivesMATLABimagesofeachsolution.Inordertofacilitatethefurtherresearchofvibrationprocess,alsogivesavibrationprocessanimation.Formorecomplexboundaryconditionsandthereforeshouldnotbeobtainedbyanalyticsolutionwaveequation,providesarelativelysimplemethod.
KEYWORDS:
Waveequation;finitedifferencemethod;showdifferencescheme;numericalsolution
目录
引言 1
1有限差分法 2
1.1 有限差分法 2
2显式差分格式 5
2.1波动方程的差分格式 5
2.2初始、边界条件的差分格式 5
3计算示例 8
3.1 边界条件为零,初始位移、初始速度不为零 8
3.2初始条件为零,边界条件不为零 12
3.3自由项不为零 17
3.4自由项、定解条件都不为零 21
总结 26
致谢 26
参考文献 26
引言
在实际工作中,经常涉及求解常微分方程边值问题、偏微分方程边值问题.通常,方程的精确解难以求解,而只能得到求解区域内某些点处解的近似值,即数值解.在数学上,二阶线性偏微分方程分为双曲型、抛物型、椭圆型三种,其中双曲型方程的典型代表是波动方程.当满足合适的边界条件和初始条件时,波动方程是可以求解的.常用的求解方法有分离变量法[1,2,3,5,6,7]、行波法[1,3,5,6]、格林函数法[1,3,4]、积分变换法[1,2,3]等等.其中最为常用的是分离变量法.而当所讨论的问题的边界条件较为复杂时,分离变量法将不再便于计算.因此,对于不易求解解析解的波动方程的数值解的研究和发展具有重要意义.
有限差分法是一种较为成熟的数值解法,它的主要思想是将连续的问题离散化.是将每一处导数由有限差分近似公式代替,从而把求解偏微分方程的问题转换成求解代数方程的问题的方法.本文在介绍了有限差分法的基础知识后给出了波动方程及其初始、边界条件的差分格式,列举了一些不同定解条件下的波动方程的数值求解方法及图像.
1有限差分法
1.1有限差分法
有限差分法是在采用数值计算方法求解偏微分方程时,将每一处导数由有限差分近似公式代替,从而把求解偏微分方程的问题转换成求解代数方程的问题的方法.有限差分法首先用差分代替微分方程中的微分,将连续变化的变量离散化,从而得到差分方程组的数学形式;继而再求解差分方程组.
设函数为,其独立变量有一很小的增量,则该函数的增量为:
(1.1)
称为函数的一阶差分,它与微分不同,因是有限量的差,故称为有限差分.而一阶差分除以增量的商,即为一阶差商:
(1.2)
一阶差分仍是独立变量的函数,类似地,可以得到二阶差分,即
(1.3)
图2.1.1差分运算原理图
显然,只要上述增量很小,差分与微分之间的差异将很小.由于一阶导数
(1.4)
是无限小的微分除以无限小微分商.因此,如图2.1.1所示,当用2、3点的量来表示函数的导数时,可近似地用差分表示为:
(向前差分)(1.5)
即用有限差分除以增量的商来表示,也称为差商.同理,当用1、2点的量来表示函数的一阶导数时,可近似表达为:
(向后差分)(1.6)
而用1、3点的量来表示时,则有:
(中心差分)(1.7)
可见,函数在处的导数可用三种差商近似表示,但哪一种精确度更高,则可以通过泰勒公式的展开式来分析.
(1.8)
和
(1.9)
可见,将(1.8)式中的二阶以及二阶以上项略去后,便可得到,这就是(1.5)式.同理,由(1.9)式可得到(1.6)式,它们都截断于项,而把项以及更高幂次的项全部略去.(1.7)式相当于把相应的泰勒公式(1.8)式及(1.9)式联立起来,得到:
(1.10)
当略去项以及更高幂次的项以后,便得到
即(1.7)式.很明显,以上三种差商表达式中,以(1.7)式所示的差商的截断误差最小,因为(1.7)式略去的是三阶以上的无穷小,而其他的是二阶以上的无穷小.
对于二阶导数,同样可用一阶差商的差商来近似表示,即
(1.11)
上式相当于把泰勒公式
(1.12)
截断到项,略去了项以及更高幂次的项.
偏导数也可以仿照上述方法表示为差商,它用各离散点上函数的差商来近似代替该点的偏导数,将需求解的边值问题转化为一组差分方程,而后根据差分方程组(代数方程组),解出位于各离散点上的待求函数值,便可得到相应的数值解.
由上可见,有限差分法是以差分原理为基础的一种数值方法,它实质上是将物理上连续域的问题变换为离散系统的问题来求解,也就是用网格状离散化模型上各离散点的数值解来逼近连续场域的真实解的.
有限差分法的应用范围很广,不但能求解均匀媒质和不均匀线性媒质中的位场,而且还能求解非线性媒质中的场;它不仅能求解恒定场或似稳场,还能求解时变场.在边值问题的数值方法中,此法是相当简便的,在计算机存储容量容许的情况下,有可能采用较精细的网格,使离散化模型能较精细地逼近真实问题,获得具有足够精度的数值解.
2显式差分格式
2.1波动方程的差分格式
考虑如下一维波动方程的定解问题
(2.1)
取进行离散化,节点坐标为
(2.2)
节点处的函数为.在点,、用中心差商代替,则式(1.1)中的偏微分方程变为:
(2.3)
引入,上式变为
(2.4)
由和时刻的可直接求时刻的,不必解联立方程组,故这种差分格式是显式的.误差为,可以证明,当时,这种格式是收敛和稳定的.
2.2初始、边界条件的差分格式
初始条件用向前差商代替,有
(2.5)
初始条件变为
(2.6)
为找到误差为边界条件的差分式,取等分坐标轴,结点坐标为,相应的有,则有泰勒展开
(2.7)
(2.8)
(2.9)
(2.10)
由(2.7)可得一阶向前差商公式
(2.11)
由(2.8)可得一阶向后差商公式
(2.12)
式(2.9)式(2.7),可得二阶向前差商公式
(2.13)
式(2.10)式(2.8),可得二阶向后差商公式
(2.14)
将二阶向前差商公式(2.13)代入(2.7),则:
由此可得一阶向前差商
(2.15)
同理将二阶向前差商公式(2.14)代入(2.8),可得一阶向后差商
(2.16)
将(2.15)代入边界条件后得
令得到左边界条件为:
(2.17)
将(2.16)代入边界条件同理可得
令得到右边界条件为:
(2.18)
3计算示例
3.1边界条件为零,初始位移、初始速度不为零
计算一维波动方程
微分方程的差分格式为
其中令,,因此有(其中单位为,单位为,单位为,单位为),上式变为
其节点坐标为
其中,并选.
其初始条件
可写为
边界条件
由于,则边界条件可写为
计算程序如下:
%计算程序
clear
%设置参数
c=1;
P=1;
L=1;
h=0.05;
T=1;
tao=P*h/c;
N=fix(L/h);
图3.1此振动图像是时间每间隔为0.05秒时的位移随坐标的变化图像.
M=fix(T/tao);
%设置u矩阵及x的值
I=N+1;
K=M+1;
forii=1:
I
x(ii)=(ii-1)*h;
end
u(I,K)=zeros;
%设置初始条件
forii=1:
I
u(ii,1)=sin((ii-1).*h.*pi);
u(ii,2)=sin((ii-1).*h.*pi)+(ii-1).*h.*tao.*(1-(ii-1).*h);
end
%设置左端第一类边界条件
u(1,:
)=0;
%设置右端第一类边界条件
u(I,:
)=0;
%计算矩阵u
fork=2:
K-1
forii=2:
I-1
u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1);
end
end
u;
fork=1:
K
plot(x,u(:
k),'-k','LineWidth',2)
holdon
end
%axis([0,1,0,1])
xlabel('\fontsize{14}\bfx/cm')
ylabel('\fontsize{14}\bfu/cm')
gridon
图3.2此图像为图3.1所表示的振动过程中初始时刻(用点划线表示)及时间为0.25s时刻(用实线表示)的波动图像.
图3.3此图像为图3.1所表示的振动过程中0.25s时刻(用点划线表示)及时间为0.5s时刻(用实线表示)的波动图像.
图3.4此图像为图3.1所表示的振动过程中0.5s时刻(用点划线表示)及时间为0.75s时刻(用实线表示)的波动图像.
图3.5此图像为图3.1所表示的振动过程中0.75s时刻(用点划线表示)及时间为1s时刻(用实线表示)的波动图像.
计算结果如图3.1所示,图中曲线至上而下时间间隔为0.05秒.将计算程序稍作修改,物理过程图像如下:
此过程可以看做是两端固定初始时刻位移为,初始速度为的弦的振动,具体过程分析如下:
、如图3.2所示,点划线表示振动过程的初始时刻的状态,实线表示振动过程的时刻的状态.初始时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
、如图3.3所示,点划线表示振动过程的时刻的状态,实线表示振动过程的时刻的状态.时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
、如图3.4所示,点划线表示振动过程的时刻的状态,实线表示振动过程的时刻的状态.时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
、如图3.5所示,点划线表示振动过程的时刻的状态,实线表示振动过程的时刻的状态.时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
为观察实际振动时的情况,可将上面程序修改为如下动画程序:
clear
图3.6此图像为振动过程中的截图.
%设置参数
c=1;
P=1;
L=1;
h=0.05;
T=1;
tao=P*h/c;
N=fix(L/h);
M=fix(T/tao);
%设置u矩阵及x的值
I=N+1;
K=M+1;
forii=1:
I
x(ii)=(ii-1)*h;
end
u(I,K)=zeros;
%设置初始条件
forii=1:
I
u(ii,1)=sin((ii-1).*h.*pi);
u(ii,2)=sin((ii-1).*h.*pi)+(ii-1).*h.*tao.*(1-(ii-1).*h);
end
%设置左端第一类边界条件
u(1,:
)=0;
%设置右端第一类边界条件
u(I,:
)=0;
%计算矩阵u
fork=2:
K-1
forii=2:
I-1
u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1);
end
end
u;
fork=1:
K
plot(x,u(:
k),'-k','LineWidth',2)
holdon
comet(x,u(:
k))
pause(0.5)
end
%axis([0,1,0,1])
xlabel('\fontsize{14}\bfx/cm')
ylabel('\fontsize{14}\bfu/cm')
gridon
3.2初始条件为零,边界条件不为零
计算一维波动方程:
令(其中单位为,单位为,单位为,单位为),编写计算程序如下.
%计算程序
clear
%设置参数
c=1;
P=1;
L=1;
h=0.05;
T=1;
tao=P.*h/c;
N=fix(L/h);
图3.7此振动图像是时间每间隔为0.05秒时的位移随坐标的变化图像.起始时为一横轴上的一条直线,最后为图像中最左边的曲线.
M=fix(T/tao);
%设置u矩阵及x的值
I=N+1;
K=M+1;
forii=1:
I
x(ii)=(ii-1)*h;
end
u(I,K)=zeros;
%设置初始条件
forii=1:
I
u(ii,1)=0;
u(ii,2)=0;
end
%设置左端第一类边界条件
u(1,:
)=0;
%设置右端第一类边界条件
forii=1:
K
u(I,ii)=1.*sin(2.*pi.*(ii-1).*tao);
end
%计算矩阵u
fork=2:
K-1
forii=2:
I-1
u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1);
end
end
u;
fork=1:
K
plot(x,u(:
k),'-k','LineWidth',2)
holdon
end
%axis([0,1,0,1])
xlabel('\fontsize{14}\bfx/cm')
ylabel('\fontsize{14}\bfu/cm')
gridon
图3.8此图像为图3.7所表示的振动过程中初始时刻(用点划线表示)及时间为0.25s时刻(用实线表示)的波动图像.
图3.9此图像为图3.7所表示的振动过程中0.25s时刻(用点划线表示)及时间为0.5s时刻(用实线表示)的波动图像.
图3.10此图像为图3.7所表示的振动过程中0.5s时刻(用点划线表示)及时间为0.75s时刻(用实线表示)的波动图像.
图3.11此图像为图3.7所表示的振动过程中0.75s时刻(用点划线表示)及时间为1s时刻(用实线表示)的波动图像.
计算结果如图3.7所示,起始时间曲线的位置与横轴吻合,最终的振动位移如图中最左面的曲线.其中为间隔0.05秒时的振动曲线.将计算程序稍作修改,物理过程图像如下:
此振动过程可以看做是始时刻的位移为零,初始速度为零,左端固定,右端位移为的弦的振动,具体过程分析如下:
、如图3.8所示,点划线表示振动过程的初始时刻的状态,实线表示振动过程的时刻的状态.初始时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
、如图3.9所示,点划线表示振动过程的时刻的状态,实线表示振动过程的时刻的状态.时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
、如图3.10所示,点划线表示振动过程的时刻的状态,实线表示振动过程的时刻的状态.时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
、如图3.11所示,点划线表示振动过程的时刻的状态,实线表示振动过程的时刻的状态.时刻各个质点波动图像为点划线形状,经过后,各个质点波动图像变为实线形状.
为观察实际振动时的情况,可将上面程序修改为如下动画程序:
clear
%设置参数
v=1;
图3.12此图为振动过程中的截图.
P=1;
L=1;
h=0.02;
T=1;
tao=P*h/v;
N=fix(L/h);
M=fix(T/tao);
%设置u矩阵及x的值
I=N+1;
K=M+1;
forii=1:
I
x(ii)=(ii-1)*h;
end
u(I,K)=zeros;
%设置初始条件
forii=1:
I
u(ii,1)=0;
u(ii,2)=0;
end
%设置左端第一类边界条件
u(1,:
)=0;
%设置右端第一类边界条件
forii=1:
K
u(I,ii)=1.*sin(2*pi.*(ii-1).*tao);
end
%计算矩阵u
fork=2:
K-1
forii=2:
I-1
u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1);
end
end
u;
fork=1:
K
axis([0,1,-1,1])
plot(x,u(:
k),'-k','LineWidth',2)
holdon
comet(x,u(:
k))
pause(0.5)
end
xlabel('\fontsize{14}\bfx/cm')
ylabel('\fontsize{14}\bfu/cm')
gridon
3.3自由项不为零,初始、边界条件为零
计算一维波动方程
令(其中单位为,单位为,单位为,单位为),编写计算程序如下.
%计算程序
clear
%设置参数
图3.13此振动图像是时间每间隔为0.05秒时的位移随坐标的变化图像.起始时为一横轴上的一条直线,最后为图像中左部最下面的曲线,右部为最上面的曲线.
A=1;
omega=2*pi;
c=1;
P=1;
L=1;
h=0.05;
T=1;
tao=P*h/c;
N=fix(L/h);
M=fix(T/tao);
%设置u矩阵及x的值
I=N+1;
K=M+1;
forii=1:
I
x(ii)=(ii-1)*h;
end
u(I,K)=zeros;
%设置初始条件
forii=1:
I
u(ii,1)=0;
u(ii,2)=0;
end
%设置左端第一类边界条件
u(1,:
)=0;
%设置右端第一类边界条件
u(I,:
)=0;
%计算矩阵u
fork=2:
K-1
forii=2:
I-1
u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1)...
+tao^2.*A.*cos((ii-1)*h.*pi/L).*sin(omega.*(k-1).*tao);
end
end
u;
fork=1:
K
plot(x,u(:
k),'-k','LineWidth',2)
holdon
end
xlabel('\fontsize{14}\bfx/cm')
ylabel('\fontsize{14}\bfu/cm')
gridon
计算结果如图3.13所示,起始时为一横轴上的一条直线,最后为图像中左部最下面的曲线,右部为最上面的曲线.
图3.14此图像为图3.13所表示的振动过程中初始时刻(用点划线表示)及时间为0.25s时刻(用实线表示)的波动图像.
图3.15此图像为图3.13所表示的振动过程中0.25s时刻(用点划线表示)及时间为0.5s时刻(用实线表示)的波动图像.
图3.16此图像为图3.13所表示的振动过程中0.5s时刻(用点划线表示)及时间为0.75s时刻(用实线表示)的波动图像.
图3.17此图像为图3.13所表示的振动过程中0.75s时刻(用点划线表示)及时间为1s时刻(用实线表示)的波动图像.
将计算程序稍
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 物理学 毕业论文 波动 方程 数值 解法 研究
![提示](https://static.bdocx.com/images/bang_tan.gif)