程冰，于加举，吴自库，陈秀荣

(青岛农业大学理学与信息科学学院，山东青岛 266109)

波动方程或波方程，是由麦克斯韦方程组导出的、描述电磁场波动特征的一组偏微分方程，主要描述自然界中的各种波动现象，例如声波、光波和水波[1-3]。最近几十年，波动方程引起了许多领域学者的广泛关注，尤其在科学工程领域，人们更加重视对波动方程解的研究。有些波动方程可以通过经典方法求得精确解，例如多级局部时间步方法等[4]。由于不同条件的约束，多数情况无法获得波动方程的精确解或解析解。与经典方法相比，数值方法有时更有效。常用的数值方法有高阶紧致差分法[1,3]、有限差分法[5]等，但这些方法多依赖于初值条件和边值条件，因此稳定性较差。

最小二乘支持向量机(least squares support vector machine，LS-SVM)方法是用训练误差的平方代替支持向量机(support vector machine，SVM)中的松弛变量，并用等式约束代替不等式约束，避免解二次规划问题，可以求得模型参数的解析解[6]。该方法已经用于求解常微分方程，取得非常好的近似效果[7]，但目前还很少用于求解偏微分方程。本文利用LS-SVM方法求解一维波动方程，该方法不依赖于边值条件和初值条件，所得近似解以标准形式给出，可以调节参数使误差最小。该近似解由两部分组成，一部分是满足边值条件的已知函数，另外一部分是两项的乘积，其中一项是在边界上取值为0的已知函数，另一项是与径向核函数相关的未知函数。

本文首先给出LS-SVM的原理和求解一维波动方程的过程，然后给出参数选择方法和稳定性分析，最后通过两个数值算例验证方法的有效性。

1 用LS-SVM方法求解一维波动方程

考虑一维波动方程：

(1)

(2)

(3)

B(V)=x(l-x)t2

(4)

将方程(2)代入方程(1)得：

(5)

其中，

Q(V)=Btt-k(x)Bxx，W(V)=Att-k(x)Bxx

G(V,Vj)=G1(V,Vj)+G2(V,Vj)+

G3(V,Vj)+G4(V,Vj)

G1(V,Vj)=Q(V)Φ(V,Vj)，G2(V,Vj)=

B(V)[Φtt-k(x)Φxx]

G3(V,Vj)=-2k(x)Bx(V)Φx，G4(V,Vj)=

2k(x)Bt(V)Φt

为求得回归参数，将样本点代入式(5)，则原问题由LS-SVM转化为求解如下二次规划问题：

(6)

约束条件：

(7)

其中，式(6)为目标函数，式(7)为约束条件，γ∈R+，为正则化参数，e是偏差项。上述优化问题的拉格朗日函数为：

(8)

依据卡罗需-库恩-塔克(Karush-Kuhn-Tucker,KKT)最优化条件有：

(9)

消去ei，回归参数由方程组(10)给出：

(10)

2 参数选择与近似解的稳定性分析

2.1 参数选择

方程组(10)含有2(M+1)个未知变量，容易解出。LS-SVM的参数对算法的性能有着非常重要的影响，参数选择不同，近似效果也不同。然而解的精度依赖于正则化因子γ及高斯核宽度参数σ。γ越大越好，一般取值不小于1×105。如何选择合适的γ和σ，目前还没有明确的理论方法，通常采用试凑法，但试凑法常使优化问题的解受主观因素的影响，且费时费力。

本文选用网络搜索法来确定模型参数，以更好地近似求解一维波动方程。首先确定γ和σ的范围，然后确定步长。以γ和σ为坐标轴，绘出误差等高线，从而确定最佳参数。

2.2 近似解的稳定性分析

(11)

显然，‖Lu(V)‖∞可以作为衡量解的稳定性和精度的指标。由于Lu(V)连续可微，且在边界上均为0，因而‖Lu(V)‖∞在内部取得。不妨令：

‖Lu(V)‖∞=|Lu(x*,t*)|,(x*,t*)∈(0,l)×(0,b)

(12)

假设距离点V(x*,t*)最近的样本点为V*(xi0,ti0)，于是有：

‖Lu(V)‖∞

=|L(x*,t*)-Lu(xi0,ti0)+Lu(xi0,ti0)|

≤|L(x*,t*)-Lu(xi0,ti0)|+|Lu(xi0,ti0)|

≤Md(τ+h)+‖e‖

(13)

这里Md为正常数，Md=max{‖Lux(V)‖∞,Lut(V)‖∞}。由(12)和(13)可知该算法是稳定的。

3 数值例子

通过两个数值例子验证该方法的有效性，并对所得近似解和精确解进行比较。式(14)定义的最大误差(Emax)和式(15)定义的平均误差(Emean)可作为衡量LS-SVM近似程度的指标。

(14)

(15)

其中uA(xi,tl)是精确解，uL(xi,tl)是由LS-SVM求得的近似解。在计算过程中，σ和γ都是常数，固定γ，研究最大误差随σ变化的过程，从而可以找到使最大误差取得最小值时的σ和γ。

例1考虑波动方程：

(16)

它的精确解u(x,t)=costsinx。在计算过程中，取A(x,t)=(1-t2)sinx，B(x,t)=x(π-x)t2，空间步长分别取h=0.157 1、0.078 5，时间步长分别取τ=0.1、0.025，取γ=1.0×108，σ=0.52。表1给出了两种情况下的数值结果，分别用D01、D02表示。图1A为方程(16)的精确解图像，图1B为方程(16)的LS-SVM近似解图像，图2为方程(16)的精确解与LS-SVM近似解之间的误差图像。可以看出，尽管所用的训练数据较少，但近似解和精确解非常接近。由于LS-SVM近似解与精确解之间的误差非常小，为了更好地呈现误差，对Emax取对数得lgEmax。图3给出了不同γ时，lgEmax随σ的变化曲线。由图3可以看出，当γ=1.0×109时，Emax相对稳定。

A.精确解；B.LS-SVM近似解。图1 例1的数值结果Fig.1 The numerical results of example 1

图2 例1的LS-SVM近似解与精确解误差Fig.2 Error δ between LS-SVM approximate solution and exact solution of example 1

图3 例1的最大误差对数Fig.3 Logarithm of maximum error of example 1

表1 近似解的两个数值结果Table 1 Two numerical results of approximate solution

例2考虑波动方程：

(17)

它的精确解u(x,t)=x+x2sinht。在计算过程中，取A(x,t)=x(1+sinht)，B(x,t)=x(1-x)t2。空间步长分别取h=0.05、0.03，时间步长分别取τ=0.10、0.05，表2给出了两种情况下的数值结果，分别用D03、D04表示。图4A为方程(17)的精确解图像，图4B为方程(17)的LS-SVM近似解图像，图5为方程(17)精确解和LS-SVM近似解之间的误差图像。可以看出，尽管所用的训练数据较少，但近似解和精确解非常接近。图6给出了不同γ时，lgEmax随σ的变化曲线。由曲线可以看出，σ<0.01时，Emax波动较大。

表2 LS-SVM近似解的两个数值结果Table 2 Two numerical results of LS-SVM approximate solution

A.精确解；B.LS-SVM近似解。图4 例2的数值结果Fig.4 The numerical results of example 2

图5 例2的LS-SVM近似解与精确解误差Fig.5 Error δ between LS-SVM approximate solution and exact solution of example 2

图6 例2的最大误差对数Fig.6 Logarithm of maximum error of example 2

4 结论

LS-SVM用于解偏微分方程，优点是对研究区域的形状、边界、网格点分布没有要求，适应性强，不依赖于边值条件和初值条件。本文用该方法求出一维波动方程的近似解，该近似解以标准形式给出，它由两部分组成，其中一部分是满足边值条件的已知函数，另一部分是两项的乘积，一项是边值为0的已知函数，另一项是与核函数相关的未知函数。为使误差最小，采用网络搜索法确定正则化因子和核宽度参数的取值，并证明了近似解的稳定性。两个数值例子表明，该方法求解一维波动方程具有较好的稳定性和较高的精度。因此该方法可用于求解边界更为复杂的波动方程。