基于五点中心差分算法求波动方程的解

张 琪

(山西职业技术学院 基础教学部，太原 030006)

波动方程是重要的一类偏微分方程，它是自然界中物质波动现象的数量刻画，如声波、水波、光波、电磁波、冲击波等[1]。波动方程的求解大致分2类，其一为求方程的精确解，揭示波传播中某点的性质与运动状态。郑明等[2]采用第二类Chebyshev小波方法对波动方程进行数值求解并获得较好的精度，为Haar小波的数值求解提供解题思路。但是由于测量存在误差，导致波动方程的精确解不能完整地反映实际问题[3]。因此，波动方程的研究转向了另一个方向，即求方程的近似解(数值解)，只要方程的近似解满足实际问题的精度要求，数值解就赋予了现实的意义。许小勇等[4]利用移位的第一类Chebyshev多项式建立求解费巨擘守恒条件下的波动方程数值解的方法，提高方程组转化求解的精度。提高波动方程数值求解的精度是方程求解的目标，在数值求解方法中有三点中心差分法、谱方法、有限元法等[5]。本文在经典的三点中心差分基础上选用五点中心差分法，提高数值求解的精度，精度对于位置变量由二阶提高到四阶，该算法格式清晰，易迭代、易推广、易理解。

1 三点中心差分法与五点中心差分法概述

1.1 三点中心差分法基本思想

三点中心差分法的基本思想是将求解区域网格化，在网格节点处，用有限个网格节点的函数值代替求解区域内解函数的值，然后在网格节点上用差分方程的解近似代替微分方程的解，直接求解得出基本方程和相应的定解条件的近似解。

1)三点中心差分方程原理

建立差分方程原理：利用有限差分法离散波动方程，方法多种多样而且对于同一方程可以建立不同的差分方法，同一方法可用于不同的方程。

设波动方程为：

(1)

式中：0≤x，y≤1；a为t时刻的传播速度；u(x，y，t)为在t时刻空间位置(x，y)点的位移；τ为时间步长；h为空间步长。网格点记为(xj，yk，tn)，其中xj=jh，yk=kh，tn=nτ，k，j=0，1，…，N，h=1/N，n≥0。

在建立之前要用到泰勒展开式，f(x)在x0处展开时的表达式为：

(2)

2)三点中心差分算法推导过程

导数的差分近似可以通过式(2)来计算。

(3)

(4)

由式(2)和式(3)得：

(5)

(6)

(7)

(8)

由式(7)和式(8)得：

(9)

(10)

(11)

(12)

由式(11)和式(12)得：

(13)

(14)

(15)

则结合式(14)和式(15)得：

(16)

综合式(6)和式(10)、式(16)得

(17)

取时间步长Δt=τ，空间步长Δx=Δy=h，则式(17)变为：

(18)

式(18)即为三点中心差分方程的显格式。

1.2 五点中心差分算法的推导过程

(19)

式中c1，c2，c3，c4，c5待定，将式(19)右端各项在x=jΔx处泰勒展开得：

ο(Δx6)。

(20)

ο(Δy6)。

(21)

(22)

(23)

由式(22)和式(23)得：

(24)

(25)

结合式(15)和式(25)得：

(26)

结合式(20)、式(21)和式(26)得：

(27)

取时间步长Δt=τ，空间步长Δx=Δy=h，则式(27)变为：

(28)

式(28)则为五点中心差分方程的显示格式。

2 五点中心差分算法分析

2.1 初始条件离散

因为该显示差分格式的时间层是3层的，换句话说每个时间步推进都需要知道前2个时间步的值[12]。

下面计算第1个时间步的值：

在式(28)中令n=0，则有

(29)

(30)

(31)

ν(xj，yk)τ。

(32)

2.2 波动方程数值算法步骤

2.3 五点中心差分与三点中心差分算法比较分析

在求解双曲型微分方程时，通常先将方程离散化，在每一分割点处用邻近点的函数值替代函数的导数，把求解偏微分方程的问题转换成代数方程问题，即通过有限差分法进行方程求解。求解过程差分节点数选取越多，近似解的精度越高。

1)三点中心差分法显示格式：

2)五点中心差分法显示格式：

3 结论

本文采用泰勒展开的方法构造五点中心差分法算法格式并得出初始条件的离散，目的在于提高波动方程的数值解精度，精度由位置变量的二阶提高到四阶，时间变量的精度保持不变。该算法具有格式简单、计算速度快等优点，是有一种有效且有广泛应用前景的方法。