物理学毕业论文波动方程的数值解法研究

1、 学 院 本科毕业论文（设计）题 目 波动方程的数值解法研究 院 系 物理与电子工程学院 专 业 物理学 姓 名 学 号 学习年限 2007年9月至2011年7月 指导教师 申请学位 理学 学士学位 二零一一年五月二十三日波动方程的数值解法研究摘 要: 本文利用有限差分法及其显示差分格式表示了波动方程的差分格式以及初始、边界条件的差分格式, 进而讨论了不同初始、边界条件下的波动方程的解的情况, 并利用MATLAB数学软件编程绘制出了各种解的图像. 为了方便对振动过程的进一步研究, 还给出了振动过程的动画演示. 对于方程初始、边界条件比较复杂因而无法求得解析解的波动方程, 有限差分方法的数值解法

2、给出了一种较为普适而精确的计算方法.关键词: 波动方程; 有限差分法; 显式差分格式; 数值解Numerical Solution of Wave EquationABSTRACT: In this paper, using finite difference method and show difference scheme showed that the difference scheme for wave equation and the initial and boundary conditions, then discussed the different initial and b

3、oundary conditions of the wave equation under the circumstances, and using mathematical software program gives MATLAB images of each solution. In order to facilitate the further research of vibration process, also gives a vibration process animation. For more complex boundary conditions and therefor

4、e should not be obtained by analytic solution wave equation, provides a relatively simple method.KEYWORDS：Wave equation; finite difference method; show difference scheme; numerical solution目 录引 言11 有限差分法21.1有限差分法22 显式差分格式52.1 波动方程的差分格式52.2 初始、边界条件的差分格式53 计算示例83.1边界条件为零, 初始位移、初始速度不为零83.2 初始条件为零, 边界条件

5、不为零123.3 自由项不为零173.4 自由项、定解条件都不为零21总 结26致 谢26参考文献26 引 言在实际工作中,经常涉及求解常微分方程边值问题、偏微分方程边值问题. 通常,方程的精确解难以求解,而只能得到求解区域内某些点处解的近似值, 即数值解. 在数学上, 二阶线性偏微分方程分为双曲型、抛物型、椭圆型三种, 其中双曲型方程的典型代表是波动方程. 当满足合适的边界条件和初始条件时, 波动方程是可以求解的. 常用的求解方法有分离变量法1,2,3,5,6,7、行波法1,3,5,6、格林函数法1,3,4、积分变换法1,2,3等等. 其中最为常用的是分离变量法. 而当所讨论的问题的边界条件

6、较为复杂时, 分离变量法将不再便于计算. 因此, 对于不易求解解析解的波动方程的数值解的研究和发展具有重要意义.有限差分法是一种较为成熟的数值解法, 它的主要思想是将连续的问题离散化. 是将每一处导数由有限差分近似公式代替, 从而把求解偏微分方程的问题转换成求解代数方程的问题的方法. 本文在介绍了有限差分法的基础知识后给出了波动方程及其初始、边界条件的差分格式,列举了一些不同定解条件下的波动方程的数值求解方法及图像.1 有限差分法1.1 有限差分法有限差分法是在采用数值计算方法求解偏微分方程时, 将每一处导数由有限差分近似公式代替, 从而把求解偏微分方程的问题转换成求解代数方程的问题的方法.

7、有限差分法首先用差分代替微分方程中的微分, 将连续变化的变量离散化, 从而得到差分方程组的数学形式; 继而再求解差分方程组.设函数为, 其独立变量有一很小的增量, 则该函数的增量为: (1.1)称为函数的一阶差分, 它与微分不同, 因是有限量的差, 故称为有限差分. 而一阶差分除以增量的商, 即为一阶差商: (1.2)一阶差分仍是独立变量的函数, 类似地,可以得到二阶差分, 即 (1.3)图2.1.1 差分运算原理图显然, 只要上述增量很小, 差分与微分之间的差异将很小. 由于一阶导数 (1.4)是无限小的微分除以无限小微分商. 因此, 如图所示, 当用2、3点的量来表示函数的导数时, 可近似

8、地用差分表示为: (向前差分) (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)式及

9、(1.9)式联立起来, 得到: (1.10)当略去项以及更高幂次的项以后, 便得到即(1.7)式. 很明显, 以上三种差商表达式中, 以(1.7)式所示的差商的截断误差最小, 因为(1.7)式略去的是三阶以上的无穷小, 而其他的是二阶以上的无穷小.对于二阶导数, 同样可用一阶差商的差商来近似表示, 即 (1.11)上式相当于把泰勒公式 (1.12)截断到项, 略去了项以及更高幂次的项.偏导数也可以仿照上述方法表示为差商, 它用各离散点上函数的差商来近似代替该点的偏导数, 将需求解的边值问题转化为一组差分方程, 而后根据差分方程组(代数方程组), 解出位于各离散点上的待求函数值, 便可得到相应的

10、数值解.由上可见, 有限差分法是以差分原理为基础的一种数值方法, 它实质上是将物理上连续域的问题变换为离散系统的问题来求解, 也就是用网格状离散化模型上各离散点的数值解来逼近连续场域的真实解的.有限差分法的应用范围很广, 不但能求解均匀媒质和不均匀线性媒质中的位场, 而且还能求解非线性媒质中的场; 它不仅能求解恒定场或似稳场, 还能求解时变场. 在边值问题的数值方法中, 此法是相当简便的, 在计算机存储容量容许的情况下, 有可能采用较精细的网格, 使离散化模型能较精细地逼近真实问题, 获得具有足够精度的数值解.2 显式差分格式2.1 波动方程的差分格式考虑如下一维波动方程的定解问题 (2.1)

11、取进行离散化, 节点坐标为 (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)可得一阶向后差商公式

12、 (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 边界条件为零, 初始位移、初始速度不为零计算一维波动方程微分方程的差分格式为其中令, , 因此有 (其中单位为,单位为,单位为,单位为

13、), 上式变为其节点坐标为其中, 并选.其初始条件可写为边界条件由于, 则边界条件可写为计算程序如下:%计算程序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;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%设置初始条件for ii=1:Iu(ii,1)=sin(ii-1).*h.*pi);u(ii,2)=sin(ii-1).*h.*pi)+(ii-1).*

14、h.*tao.*(1-(ii-1).*h);end%设置左端第一类边界条件u(1,:)=0;%设置右端第一类边界条件u(I,:)=0;%计算矩阵ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); endendu;for k=1:K plot(x,u(:,k),'-k','LineWidth',2) hold onend%axis(0,1,0,1)xlabel ('fontsize14bfx/cm')ylabel ('fontsize14bfu/cm'

15、)grid on图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秒. 将计算程序稍作修改, 物理过程图像如下:此过程

16、可以看做是两端固定初始时刻位移为, 初始速度为的弦的振动,具体过程分析如下:、如图3.2所示, 点划线表示振动过程的初始时刻的状态, 实线表示振动过程的时刻的状态.初始时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.3所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.4所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.5所示,点划线表示振动

17、过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.为观察实际振动时的情况, 可将上面程序修改为如下动画程序: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;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%设置初始条件for ii=1:Iu(ii,1)=sin(ii-1).*h.*pi);u(ii,2)=sin(ii

18、-1).*h.*pi)+(ii-1).*h.*tao.*(1-(ii-1).*h);end%设置左端第一类边界条件u(1,:)=0;%设置右端第一类边界条件u(I,:)=0;%计算矩阵ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); endendu;for k=1: K plot(x,u(:,k),'-k','LineWidth',2) hold on comet(x,u(:,k) pause(0.5)end%axis(0,1,0,1)xlabel ('fontsize1

19、4bfx/cm')ylabel ('fontsize14bfu/cm')grid on3.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;for ii=1:I x(ii)=(ii-1)*h;

20、endu(I,K)=zeros;%设置初始条件for ii=1:Iu(ii,1)=0;u(ii,2)=0;end%设置左端第一类边界条件u(1,:)=0;%设置右端第一类边界条件for ii=1:Ku(I,ii)=1.*sin(2.*pi.*(ii-1).*tao);end%计算矩阵ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); endendu;for k=1:K plot(x,u(:,k),'-k','LineWidth',2) hold onend%axis(0,1,0,1

21、)xlabel ('fontsize14bfx/cm')ylabel ('fontsize14bfu/cm')grid on图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时刻(用点划线表示)及时间为1

22、s时刻(用实线表示)的波动图像.计算结果如图3.7所示, 起始时间曲线的位置与横轴吻合, 最终的振动位移如图中最左面的曲线. 其中为间隔0.05秒时的振动曲线. 将计算程序稍作修改, 物理过程图像如下:此振动过程可以看做是始时刻的位移为零, 初始速度为零, 左端固定, 右端位移为的弦的振动, 具体过程分析如下:、如图3.8所示, 点划线表示振动过程的初始时刻的状态, 实线表示振动过程的时刻的状态. 初始时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.9所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状

23、, 经过后, 各个质点波动图像变为实线形状.、如图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

24、);%设置u矩阵及x的值I=N+1;K=M+1;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%设置初始条件for ii=1:Iu(ii,1)=0;u(ii,2)=0;end%设置左端第一类边界条件u(1,:)=0;%设置右端第一类边界条件for ii=1:Ku(I,ii)=1.*sin(2*pi.*(ii-1).*tao);end%计算矩阵ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); endendu;for k=1:K axis(0,1,-1,1) plot(x,u(

25、:,k),'-k','LineWidth',2) hold on comet(x,u(:,k) pause(0.5)endxlabel ('fontsize14bfx/cm')ylabel ('fontsize14bfu/cm')grid on3.3 自由项不为零, 初始、边界条件为零计算一维波动方程令 (其中单位为,单位为,单位为,单位为), 编写计算程序如下.%计算程序clear%设置参数图3.13 此振动图像是时间每间隔为0.05秒时的位移随坐标的变化图像.起始时为一横轴上的一条直线, 最后为图像中左部最下面的曲线,右部为最

26、上面的曲线.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;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%设置初始条件for ii=1:Iu(ii,1)=0;u(ii,2)=0;end%设置左端第一类边界条件u(1,:)=0;%设置右端第一类边界条件u(I,:)=0;%计算矩阵ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1) .+ta

27、o2.*A.*cos(ii-1)*h.*pi/L).*sin(omega.*(k-1).*tao); endendu;for k=1: K plot(x,u(:,k),'-k','LineWidth',2) hold onendxlabel ('fontsize14bfx/cm')ylabel ('fontsize14bfu/cm')grid on计算结果如图3.13所示, 起始时为一横轴上的一条直线, 最后为图像中左部最下面的曲线, 右部为最上面的曲线.图3.14 此图像为图3.13所表示的振动过程中初始时刻(用点划线表示)及时

28、间为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时刻(用实线表示)的波动图像.将计算程序稍作修改, 物理过程图像如下:此振动过程可以看做是两端固定, 初始时刻的位移为零, 初始速度也为零的弦在驱动力的作用下的振动, 具体过程分析如下:、如图3.14所示, 点划线表示振

29、动过程的初始时刻的状态, 实线表示振动过程的时刻的状态. 初始时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.15所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.16所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.17所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后,

30、各个质点波动图像变为实线形状.为观察实际振动时的情况, 可将上面程序修改为如下动画程序:%计算3动画程序clear 图3.18 此图为振动过程中的截图.%设置参数A=1;omega=2*pi;c=1;P=1;L=1;h=0.02; T=1;tao=P*h/c;N=fix(L/h);M=fix(T/tao);%设置u矩阵及x的值I=N+1;K=M+1;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%设置初始条件for ii=1:Iu(ii,1)=0;u(ii,2)=0;end%设置左端第一类边界条件u(1,:)=0;%设置右端第一类边界条件u(I,:)=0;

31、%计算矩阵ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1)+. tao2.*A.*cos(ii-1)*h.*pi/L).*sin(omega.*(k-1).*tao); endendu;for k=1:K axis(0,1,-0.08,0.08) plot(x,u(:,k),'-k','LineWidth',2) hold on comet(x,u(:,k) pause(0.1)endxlabel ('fontsize14bfx/cm')ylabel ('f

32、ontsize14bfu/cm')grid on3.4 自由项、定解条件都不为零考虑如下波动方程：令(其中单位为,单位为,单位为,单位为), 则编写计算程序如下:%计算程序clear%设置参数图3.19 此振动图像是时间每间隔0.05秒位移随坐标的变化图像.初始时刻为横轴上一条直线,终止时为图像左下方的曲线.c=1;P=1;L=1;h=0.05;T=1;omega=2*pi;tao =P*h/c;N=fix(L/h);M=fix(T/tao);%设置u矩阵及x的值I=N+1;K=M+1;for ii=1:Ix(ii)=(ii-1).*h;endu(I,K)=zeros;%设置初始条件f

33、or ii=1:Iu(ii,1)=0;u(ii,2)=(ii-1).*h.*(1-(ii-1).*h);endfor ii=1:K%设置左端第一类边界条件u(1,ii)=0;%设置右端第一类边界条件u(I,ii)=sin(omega.*(ii-1).*tao);end%计算矩阵ufor k=2:K-1 for ii=2:I-1u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1)+ tao2.*(ii-1).*h/(1+(ii-1).*h)2)2;endendu;for k=1:Kplot(x,u(:,k),'-k','linewidth'

34、;,2)hold onendxlabel('fontsize14bfx/cm') ylabel('fontsize14bfu/cm') grid on图3.20此图像为图3.19所表示的振动过程中初始时刻(用点划线表示)及时间为0.25s时刻(用实线表示)的波动图像.图3.21 此图像为图3.19所表示的振动过程中0.25s时刻(用点划线表示)及时间为0.5s时刻(用实线表示)的波动图像.图3.22 此图像为图3.19所表示的振动过程中0.5s时刻(用点划线表示)及时间为0.75s时刻(用实线表示)的波动图像.图3.23 此图像为图3.19所表示的振动过程中0.

35、75s时刻(用点划线表示)及时间为1s时刻(用实线表示)的波动图像.将计算程序稍作修改, 物理过程图像如下:此振动过程可以看做是左端固定, 右端位移为, 初始时刻的位移为零, 初始速度为的弦在驱动力的作用下的振动, 具体过程分析如下:、如图3.20所示, 点划线表示振动过程的初始时刻的状态, 实线表示振动过程的时刻的状态. 初始时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.21所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.22所示, 点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.、如图3.23所示,点划线表示振动过程的时刻的状态, 实线表示振动过程的时刻的状态. 时刻各个质点波动图像为点划线形状, 经过后, 各个质点波动图像变为实线形状.为观察实际振动时的情况, 可将上面