DOI: 10.11817/j.issn.1672-7207.2015.12.026

Daubechies小波有限元联系系数计算的广义最小二乘-复化梯形求积法

冯德山，王珣，杨炳坤，杜华坤

(中南大学 地球科学与信息物理学院，有色资源与地质灾害探查湖南省重点实验室，湖南 长沙，410083)

摘要：针对目前Daubechies小波有限元联系系数计算中附加方程复杂、计算结果会随附加方程个数的不同而改变、计算难度大，计算精度不高等现状，基于二维可分离小波理论，将复化梯形求积法与传统方法中未加载附加条件的滤波系数方程组相结合，先用复化梯形求积法求出若干个联系系数的初值，再结合广义最小二乘法理论，推导出基于[0,1]区间任意尺度下的Daubechies小波有限元联系系数的计算公式。研究结果表明：所提出的广义最小二乘-复化梯形求积法(GLS-CTQM)降低了联系系数求解难度，不仅能求得高精度的联系系数，而且可根据实际精度需求灵活地改变求积步长、选取不同求积公式，易于编程实现，计算效率和计算精度较梯形求积方法都有所   提高。

关键词：Daubechies小波；小波有限元；联系系数；最小二乘法；复化梯形求积法

中图分类号：O241.82；O242.21             文献标志码：A         文章编号：1672-7207(2015)12-4578-06

Generalized least squares - compound trapezoid quadrature method applied in connection coefficient computation of daubechies wavelet finite element method

FENG Deshan,WANG Xun, YANG Bingkun, DU Huakun

(Key Laboratory of Non-ferrous Resources and Geological Detection of Hunan Province,

School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

Abstract: Considering that Daubechies wavelet finite element connection coefficient is calculated by adding complex equations currently, calculation results vary with different numbers of additional equations with high difficulty, and calculation results are not precise enough, the compound trapezoid quadrature method (CTQM) was combined with filter coefficient equations without additional conditions in traditional method based on two-dimensional separable wavelet theory. The initial value of several connection coefficients was obtained by CTQM, and the calculation formula of Daubechies wavelet finite element connection coefficient based on interval[0,1] arbitrary scales was deduced in combination with generalized least squares (GLS) theory. The results show that GLS-CTQM reduces the difficulty of solving the connection coefficient, connection coefficient has high precision, quadrature step length can be changed flexibly, different quadrature formulas can be selected according to the actual precision demand, and the method can be easily programmed with higher computational efficiency and accuracy than trapezoid quadrature method.

Key words: Daubechies wavelet; wavelet finite element method; connection coefficient; generalized least-squares method; compound trapezoid quadrature method

有限元法自被提出以来，其理论体系日趋完善，已成为当前工程数值计算领域应用广泛的一种方法^[1]。然而，采用有限元法处理复杂的冲击波、地质裂隙等局部大梯度、强奇异性问题时，一般需要增加网格数或提升多项式阶次来确保精度，导致计算量增加，并且当数值结果达不到精度要求时，需要加密网格，原来形成的刚度矩阵不能够被网格细化后利用，浪费了计算资源。为了更好地求解此类复杂奇异性问题，一些学者将具有多尺度、多分辨特性的小波函数与有限元法相结合，产生了小波有限元^[2-3]。小波有限元法吸收了传统有限元法离散逼近的优点，同时又具有小波函数特有的多分辨特性。它将分析对象依次放入一个互相嵌套的函数空间序列中，根据实际需要由粗及细地逐步逼近信号，在不改变网格剖分的前提下提高其分辨率^[3]。当小波有限元选用的基函数为紧支撑、正交特性的小波时，所获得的系数矩阵是稀疏的，矩阵的条件数与维数无关，通过阈值运算，可以大大减少求解的计算量，从而提高算法效率。在众多小波中，Daubechies小波^[4]因为具有正交性、紧支撑性等优点，受到国内外学者的广泛关注。Daubechies 小波的尺度函数和小波函数通过两尺度方程式及小波方程构造，以数表和曲线的方式给出，均没有明确的显式数学表达式，因此，运用Daubechies小波有限元进行数值模拟时，联系系数的计算成为必不可少的核心环节。Latto等^[5]进行了无穷区间联系系数的计算，Beylkin等^[6]加深并拓展了这方面的研究；Chen等^[7]计算了[0，x]区间上的积分；Popovici^[8]基于Matlab编写了[0，x]区间上联系系数的计算程序；Monasse等^[9]通过构造新的小波方法求解了有限区间积分；Ko等^[10]研究了[0，1]区间上的积分；Chen等^[11]给出了任意尺度函数在[0，1]区间上刚度矩阵与载荷列阵的联系系数计算方法。采用上述方法计算联系系数时，由于系数矩阵不满秩，需要附加方程以得到唯一解，但需要补充的方程数目不容易确定，需要手动添加，且计算结果会随附加方程个数不同而改变，使计算过程无法完全实现程序自动化。而联系系数的求解精度及效率又直接关系到Daubechies小波有限元算法的效果。目前联系系数计算的主要问题是计算精度不高、计算难度较大、无法完全实现程序化，制约着小波有限元在数值求解偏微分方程方面的推广应用。为此，本文作者提出广义最小二乘-复化梯形求积法(generalized least squares-compound trapezoid quadrature method，简称GLS- CTQM)计算联系系数。该算法采用复化梯形公式预先求得若干个联系系数初值，并将该初值与未加载附加方程的滤波系数方程组相结合，再运用广义最小二乘法求联系系数。该求解方法理论简单易懂，可根据实际精度需求灵活地选取求积公式，易于编程实现。

1  Daubechies小波联系系数分类

20世纪90年代，Amaratunga等^[12-13]利用小波紧支撑特性和多分辨特性，将其作为试探函数，运用Galerkin法求解偏微分方程，显示了小波分析在数值计算中的巨大优越性。Daubechies小波有限元^[14]与样条小波有限元^[15]为应用最广泛的2类，其中Daubechies小波有限元因为没有显式表达式，需要求解尺度函数、尺度函数导数以及联系系数。尺度函数、尺度函数的导数求取相对容易，而对于联系系数的求解较复杂，根据求积区间及表达式的形式将联系系数归纳为5类。

Latto等^[5]定义无穷区间尺度函数及其导数乘积的积分，并将其称为联系系数。该尺度函数的内积值形式如下：

     (1)

式中：为尺度函数；d₁与d₂为导数的阶数；k与为小波平移系数；联系系数中左上标1代表第1种联系系数形式。该类联系系数的小波需要进行周期化处理，其应用也局限于无穷区间或周期边界条件的问题。然而，实际工程数值计算问题大多是基于有限域的，需将整个求解域离散为有限个单元，且每个单元域也是有限的，为此，Chen 等^[11]开展了[0，t]区间上的有限积分，其积分形式如下：

        (2)

Lin等^[16]应用半解析小波时间差分法求解Burgers方程时，定义了区间[0，2^j]上的有限积分(其中，j为小波的尺度)，其联系系数形式如下：

      (3)

Ko等^[10]研究了[0，1]区间上的积分，即第4种形式：

      (4)

Chen等^[11]在这些算法的基础上，给出了任意尺度函数在[0，1] 区间上形成刚度矩阵与载荷列阵的联系系数计算方法，即

   (5)

近年来，运用Daubechies小波有限元求解偏微分方程研究中，可将有限域的区间长度投射到[0,1]区间，故式(4)与(5)中的联系系数应用最广泛。现有的小波联系系数常限于上面几种形式，本文以第5种形式的联系系数计算为例，将式(5)中k和分别改成k₁和k₂，可得

    (6)

2  基于[0，1]区间的联系系数的求解

设x的积分区间为[0，1]，结合Daubechies小波尺度函数的支撑区间[0，2N-1]，可求得k₁和k₂的取值范围为

   (7)

为了方便公式推导，引入特征函数：

            (8)

令，则

           (9)

即

       (10)

故式(6)可化为

  (11)

Daubechies小波没有显式的数学表达式，由以下尺度函数和相应的小波函数给出^[2]：

       (12)

其中：

          (13)

N为Daubechies小波的阶数；为小波滤波系数。由式(12)可得

     (14)

将式(14)两边同求d 次导，则

 (15)

将式(15)代入式(11)，可得

 (16)

即可得如下方程组：

        (17)

式中：。由于，根据齐次方程组有非零解的条件，可知。所以，矩阵奇异，式(17)不能唯一确定。目前常用的求解方法是：引入与解空间的自由度个数相等的附加方程，从而组成恰定方程组，求解后可唯一确定准确的联系系数。对于添加的附加方程，Chen等^[11]给出如下计算式：

     (18)

用式(18)作为附加方程组时，需要求解系数和，但q和w取值范围比较广，计算结果会随附加方程个数不同而改变，难以确定哪组解是最优的。此后，陈雅琴等^[17]对附加方程数目进行了大量试算，发现需引入附加方程数目与小波消失矩及尺度函数导数阶数之和存在明显的关系，当消失矩p≤6及小波尺度V_j(j≤2)时较准确，而若计算过程中采用了高阶消失矩及较大的小波尺度，则计算误差会增大。

3  广义最小二乘-复化梯形求积法

本文提出的GLS-CTFQM算法是根据二维小波可分离理论，求得相应尺度函数、尺度函数导数值之间的乘积，根据Matlab可以很容易求得式(17)矩阵的秩，然后用所求联系系数个数减该矩阵的秩，可得到需求的联系系数的初值个数m，与文献[17]中所得结果一致。先用复化梯形积分法试求m个联系系数初值，再将这些初值添加到式(17)中，组成新的超定方程组，然后在残差平方和最小的情况下对数值解进行修正，即可求得到满足精度要求的联系系数。鉴于积分方法的相似性，本文中应用的复化梯形求积公式亦可用复化Simpson积分、Gauss-Legendre积分公式等来替代。

运用Daubechies小波有限元进行数值求解时，一般需要求解尺度函数及尺度函数导数，它们的唯一解较易求得。假设已求得Daubechies小波尺度函数及其导数，在步长足够细化的情况下，应用复化梯形求积公式理应可以求得准确的联系系数。但通过求解发现，单用复化梯形求积方法计算联系系数时，计算速度较慢，且对振荡函数积分时，存在个别联系系数计算精度较低的问题。为此，在计算式(17)的秩后，利用复化梯形求积方法与式(17)相结合，并在残差平方和最小的情况下对数值解进行修正，计算效率与精度均可满足要求。

Daubechies小波的尺度函数、尺度函数的导数值的求解可参照文献[18]中的算法进行。本文以DB3小波为例，该算法求得的尺度函数二分点处尺度函数的导数如图1所示。根据二维小波可分离理论及图1中求得的尺度函数的导数，可求得尺度函数的内积。图2所示为DB3小波尺度函数一阶导数内积。令

       (19)

式中：N为Daubechies小波阶数； 。式(6)可化为

         (20)

图1  DB3小波尺度函数一阶导数

Fig. 1  Scale function one order derivative of DB3 wavelet

图2  DB3小波尺度函数一阶导数内积

Fig. 2  Scale function one derivative inner product of DB3 wavelet

采用复化梯形公式求积算法^[19]。首先将区间[0，1]分为n个等分的子区间，在每个子区间内使用梯形公式，则

   (21)

可得

 (22)

根据Daubechies小波消失矩特性，具有N阶Daubechies小波尺度函数可以精确地表征不大于N-1阶的幂级数。在开区间内必然存在1点，使得

            (23)

故式(22)可改写为

             (24)

由于为定值，当区间[0,1]离散步长h足够小时，有。于是，式(22)可近似为

，

                 (25)

则可由式(25)试算m个平移系数k₁和k₂(一般取k₁=0，0≤k₂≤2-2N)所对应的联系系数，即

           (26)

式中：。在大多数工程实践中，，根据试算的联系系数个数与的关系，试算的附加辅助方程个数不超过8个。满足此条件时，联立式(25)与式(17)可得到如下方程组：

              (27)

式中：

E₁为将E₀扩充为的矩阵，未赋值的矩阵元素置为0，此时，对组建的新方程组(27)，在求解过程中可采用残余误差最小二乘法求得方程组的1个最近似解^[20]。设残余误差，则使即可得到方程组的一组最小二乘解，此解即为满足精度要求的联系系数。由于用复化梯形求积公式试算的联系系数计算相对简单且计算个数较少，所以，由(27)式计算联系系数时不会使计算量大增。

4  计算实例与分析

实例：的求解。

若计算实例中的联系系数，则传统Daubechies小波联系系数的求解方法一般在(17)式后附加如下方程

               (28)

即可唯一确定联系系数。若采用本文提出的GLS-CTQM方法，以步长h=1/2¹⁰为例，则可先用复合梯形积分法求出：

 (29)

再将(29)与(17)式联立，并将方程组在残差平方和最小化下求解，利用编制的Matlab程序，经0.129 s即可求得，其联系系数见图3(a)。

Landragin-Frassati等^[21]也给出了该联系系数的计算结果，若保留8位有效位数字，则这2种方法所得计算结果完全一致。只采用复化梯形求积法计算，经14.444 s可求得同样的联系系数，该联系系数及相对误差见图3(b)和3(c)。从图3(c)可见：采用复化求积法直接对振荡函数积分，在个别点存在较大误差，局部精度不高；而采用GLS-CTQM法对实例1的联系系数求解计算速度快，而且精度高。

图3  不同方法计算联系系数的结果及相对误差

Fig. 3  Values and relative error of calculating connection coefficient  by different methods

5  结论

1) 采用传统积分方法亦可求得相应的联系系数，但求解效率低，且计算结果不稳定。将复化梯形求积公式与未加载附加方程的滤波系数方程组相结合，运用最小二乘法在残差平方和最小的情况下可求得满足精度要求的联系系数。该方法求解理论易懂，可灵活地根据实际精度需求选取求积公式。

2) 本文提出的GLS-CTQM算法所得联系系数在计算精度上不但与传统的联系系数计算方法所得计算结果较吻合，而且本文提出的联系系数计算方法能有效地克服补充方程的精确数目不易确定、计算结果随附加方程个数的不同而改变、计算难度较高的问题，特别对于d₁+d₂和N较小的联系系数，计算速度非常快，更易编程实现。

参考文献：

[1] Polycarpou A C.Introduction to the finite element method in electromagnetics[M]. California:Morgan & Claypool, 2006: 1-7.

[2] 何正嘉，陈雪峰，李兵，等.小波有限元理论及其工程应用[M]. 北京: 科学出版社, 2006: 5-15.

HE Zhengjia, CHEN Xuefeng, LI Bing, et al. Theory of the wavelet based finite element methods and the application in engineering[M]. Beijing: Science Press, 2006: 5-15.

[3] 张新明, 刘克安, 刘家琦. 流体饱和多孔隙介质二维弹性波方程正演模拟的小波有限元法[J]. 地球物理学报, 2005, 48(5): 1156-1166.

ZHANG Xinming, LIU Kean, LIU Jiaqi. A wavelet finite element method for the 2-D wave equation in fluid-saturated porous media[J]. Chinese J Geophys, 2005, 48(5): 1156-1166.

[4] Daubechies I. Orthonormal basis of compactly supported wavelets[J]. Communications on Pure and Applied Mathematics, 1998, 41(XLI): 909-996.

[5] Latto A, Resnikoff L, Tenenbaum E. The evaluation of connection coefficients of compactly supported wavelets[R]. Cambridge: University of Cambridge, 1999: 1-10.

[6] Beylkin G. On the representation of operators in bases of compactly supported wavelets[J]. Communications on Pure and Applied Mathematics, 1992, 6(6): 1716-1740.

[7] Chen M Q, Hwang C, Shih Y P. The computation of Wavelet-Galerkin approximation on a bounded interval[J]. International Journal for Numerical Methods in Engineering, 1996, 39(17): 2921-2944.

[8] Popovici C I. Matlab evaluation of the coefficients for PDE solving by wavelet- Galerkin approximation[J]. Analele Stiintifice ale Universitatii Ocidius Constanta, 2010, 18(1): 287-294.

[9] Monasse P, Perrier V. Orthonormal wavelet bases adapted for partial differential equations with boundary conditions[J]. SIAM Journal on Mathematical Analysis, 1998, 29(4): 1040-1065.

[10] Ko J, Kurdila A J, Pilant M S. A class of finite element methods based on orthonormal, compactly supported wavelets[J]. Computational Mechanics, 1995, 16(4): 235-244.

[11] CHEN Xuefeng, YANG Shengjun, MA Junxing, et al. The construction of wavelet finite element and its application[J]. Finite Elements in Analysis and Design, 2004, 40(5/6): 541-554.

[12] Amaratunga K, Williams J R, Qian S, et al. Wavelet-Galerkin solutions for one-dimensional partial differential equations[J]. International Journal for Numerical Methods in Engineering, 1994, 37(16): 2703-2716.

[13] Amaratunga K, Williams J R. Wavelet based Green's function approach to 2D PDES[J]. Engineering Computations, 1993, 10(4): 349-367.

[14] Somchai S I, Eckart S. Wavelet-Galerkin solution of a partial differential equation with nonlinear viscosity[J]. Applied Mathematical Sciences, 2013, 38(7): 1849-1880.

[15] Wu C W, Chen W H. Extension of spline wavelets element method to membrane vibration analysis[J]. Computational Mechanics, 1996, 18(1): 46-54.

[16] Lin E B, Zhou X. Connection coefficients on an interval and wavelet solutions of burgers equation[J]. Journal of Computational and Applied Mathematics, 2001, 135(1): 63-78.

[17] 陈雅琴, 张宏光, 党发宁. Daubechies 小波有限元法联系系数计算研究[J]. 西安理工大学学报, 2011, 27(2): 171-176.

CHEN Yaqin, ZHANG Hongguang, DANG Faning. Research on connection coefficient computation in Daubechies wavelet finite element method[J]. Journal of Xi’an University of Technology, 2011, 27(2): 171-176.

[18] 陈雅琴. 桥梁结构分析的广义变分原理-Daubechies条件小波法研究[D]. 西安: 长安大学公路学院, 2008: 56-83.

CHEN Yaqin. Study on Generalized Variational Principle in bridge structural analysis-Daubechies conditional wavelet[D]. Xi’an: Chang’an University. School of Highway, 2008: 56-83.

[19] 李庆扬, 王能超, 易大义. 数值积分[M]. 北京: 清华大学出版社, 2009: 127-129.

LI Qingyang, WANG Nengchao, YI Dayi. Numerical analysis[M]. Beijing: Tsinghua University Press, 2009: 127-129.

[20] 魏木生. 广义最小二乘问题的理论和计算[M]. 北京: 科学出版社, 2008: 258-290.

WEI Musheng. The theory and computation of generalized least square problem[M]. Beijing: Science Press, 2008: 258-290.

[21] Landragin-Frassati A, Bonnet S, Silva A D, et al. Application of a wavelet-Galerkin method to theforward problem resolution in fluorescencediffuse optical tomography[J]. Optics Express, 2009, 17(21): 18433-18448.

(编辑  陈灿华)

收稿日期：2014-12-20；修回日期：2015-02-01

基金项目(Foundation item)：国家自然科学基金资助项目(41574116)；中南大学创新驱动项目(2015CX008)；教育部新世纪优秀人才支持计划项目(NCET-12-0551)；中南大学研究生自主探索创新项目(2015ZZTS249)；中南大学教师研究基金资助项目(2014JSJJ001)；湖湘青年创新创业平台培养对象共同资助项目(2013)(Project(41574116) supported by the National Natural Science Foundation of China；Project(2015CX008) supported by Innovation Driven of Central South University；Project(NCET-12-0551) supported by the New Century Talents of Ministry of Education；Project(2015ZZTS249) supported by Graduate Innovation and Creativity Funds of Central South University; Project(2014JSJJ001) supported by the Research Foundation of Central South University；Project(2013) supported by Hunan Youth Innovation and Entrepreneurship Training Platform)

通信作者：冯德山，教授，博士，从事地球物理正反演算法研究；E-mail：fengdeshan@126.com