用径向基函数解偏微分方程

1、用径向基函数解偏微分方程 *张颖超( 广西师范大学数学科学学院，中国 桂林541004)摘 要 讨论了用正定径向基函数解偏微分方程，通过一个数值算例，说明这个方法是可行的 针 对 数 值 算例，比较了在相同步长时，不同的正定径向基函数对微分方程数值解的精确程度，并比较不同的正定径向基函数在 相同的形状参数时绝对误差的差异，说明微分方程数值解的精确程度与径向基函数形状参数的取值密切相关 同时也论证了在插值过程中所得到的矩阵方程解的存在唯一性关键词 径向基函数; 数值解; 偏微分方程; 线性无关中图分类号 O241 81文献标识码 A文章编号 1000-2537( 2011) 05-0001-06

2、Application of Radial Basis Functions for Partial Differential EquationsZHANG Ying-chao*( School of Mathematics，Guangxi Normal University，Guilin 541004，China)AbstractAn algorithm for partial differential equations based on the positive definite radial basis funtions( RBF) approximation scheme is pre

3、sented One model problem of the algorithm is given The comparison is made with the exact solutions of the problem by different shape parameter when different radial basis functions arechosen Numerical results show that method offers a very high accuracy in computation of the partial differentialequa

4、tion It shows that choice of shape parameter is important The obtained coefficient matrix is proved to be non- singular，that is，matrix equation has a solutionKey wordsRBF; numerical solutions; partial differential equation; linear independent一般地，许多物理现象和工程技术问题都可以归结为一个微分方程 微分方程的解析解可以通过分离变量法和积分变换法等方法得到

5、 微分方程的数值逼近可以通过 Euler 法，Runge-Kutta 法和数值积分等方法得到 1-2近几十年来，人们的主要目标是寻找各种各样的无网格方法，利用径向基函数( RBF) 解微分方程是受到普遍关注的无网格方法 1971 年，Hardy3总结评论了关于 multiquadric( MQ) 函数的各种应用，特别是在地理，遥感，信号系统等方面的成功应用 自从 Kansa4-5 用径向基函数解偏微分方程( PDE) ，得到非常精 确 的解后，用径向基函数解偏微分方程引起越来越多的关注，并且 Madych 和 Nelson 3，6-7 证实了 MQ 函数插 值的收敛性近 10 年来，用径向基函

6、数配置法解偏微分方程受到广泛关注，人们已经用径向基函数配置法解线性和非线性的偏微分方程 8-11，有的用径向基函数配置法逼近椭圆型偏微分方程数值解 12-13，并且取得了不错 的结果，但是在逼近过程中所得到矩阵方程的系数矩阵是否可逆还没有被验证，即，数值解的唯一性还没有被验证* 收稿日期: 2011-02-21基金项目: 国家自然科学基金资助项目( 11171102)本文主要讨论用正定径向基函数14 解偏微分方程 ( p u) + ( p u) + qu = f( x，y) ，( x，y) ，xxyy( 1)Bu = g( x，y) ，( x，y) ，其中 p 0，q 0 或者 p 0，q 0

7、， R2 ，f( x，y) ，g( x，y) 是连续函数 本文主要讨论在插值过程中，所得到的矩阵方程的解的唯一性，还讨论了径向基函数的形状参数对数值解的影响以下是常用的正定径向基函数= e cr2 ( c 0) ;高斯径向基函数( Gauss) : ( r)1逆 MQ 函数( Inverse Multiquadrics) : ( r)( c 0) ;=槡c2 + r22c分式径向基函数( Fraction) : ( r)( c 0) =c2 + r21正定径向函数的应用1 1相关定理对于任意阶可微函数 u( r) ，r Rn ，在 内配置 N 个散乱的数据点 r ，r ，r ，令1 2NNu(

8、 r) uN ( r) = j ( r rj ) ，( 2)是正定径向基函数，· 是欧几里得范数 本j = 1其中 1 ，2 ，N 是待定系数，( r rj ) ( j= 1，2，)文取 n = 2，即 r引理 114处大于零引理 215= ( x，y) = ( x) 是正定函数的充分必要条件是 ( x) 的傅里叶变换 珦( ) 几乎处函数 ( x)若 F( x) 是 ( x) 的傅立叶变换，那么 F( x) = ( i) 2 F( x) 考虑正定径向基函数 ( r rj ) ( j = 1，2，) ，设( r rj ) = p( xx ( r rj ) + yy ( r rj )

9、) + q( r rj ) ( j= 1，2，) ，因为 r = ( x，y) ，所以 ( r rj ) = ( x xj ，y yj ) ( j = 1，2，) ，且 ( x xj ，y yj ) 的傅立叶变换是+ + = 珦( 1 ，2 ) i1( x xj) i2( y yj)( x xj ，y yj ) eedxdy 定理 1若 ( x xj ，y yj ) 的傅立叶变换是珦( 1 ，2 ) = F( x xj ，y yj ) ，那么当 p 0 且 q 0 时，珦( 1 ，2 ) 0; 当 p 0 且 q 0 时，珦( 1 ，2 )证 根据傅立叶变换的性质( 3) 0珦( 1 ，2 )

10、 = F( x xj ，y yj) = Fp( xx ( r rj ) + yy ( r rj ) ) + q( r rj ) =pFxx ( r rj ) + pFyy ( r rj ) + qF( r rj ) =22( i1 ) pF( r rj ) + ( i2 ) pF( r rj ) + qF( r rj ) =22 ( p1 + p2 q) F( r rj ) + F( r rj ) 因为 ( r rj ) ( j= 1，2，)是正定径向基函数，根据引理 1，220，珦( 1 ，2 )当 p 0 且 q 0 时， ( p1 + p2 q) q) 0; 0220，珦( ，2 )当

11、p 0 且 q 0 时， ( p1 + p21当 p 0 且 q 0，( r r1 ) ，( r r2 ) ，( r rN ) ， 线性无关定理 2NN从 ( r r1 ) ，( r r2 ) ，( r rN ) ，中取截断 ( r rN ) j = 1 假设 ( r rN ) j = 1 线性相证关，那么存在不全为零的 j ( 1 j N) 使得下式成立Nj ( r rj )j = 1= 0= ( xk ，yk ) ，那么由( 3)令 r= rk式NN+ + = j 珦( 1 ，2 ) ei1( x xj) i2( y yj)= j ( r rj )0ed1 d2= j = 1j = 1N+

12、 + ei1x ei2y 珦( ， ) e i1xj e i2yj d d=j1212 j = 1N+ + ei1xk ei2yk 珦( ， ) e i1xj e i2yj d1 d ，( 4)j122 j = 1N+ + 2 珦( ， ) d d(i1( xk xj) i2( yk yj)j e=e1212 j = 1NN+ + k ei1xk i2yk i1xj i2xj 珦( 1 ，2 ) d1 d2ej ee= k = 1j = 1NN+ + k i1xk i2yk i1xj i2yk 珦( 1 ，2 ) d1 d2 eej ee k = 1j = 1由( 4) 式，+ + N 2 珦

13、( ， ) d dj ej = 1i1( xk xj)i2( yk yj)= 0e1212 根据定理 1，N2i1( xk xj) i2( yk yj)= 0j eej = 1因此 j = 0，1 j N，所以假设不成立 所以 ( r r1 ) ，( r r2 ) ，( r rN ) ， 线性无关1 2方法分析考虑偏微分方程 ( p u) + ( p u) + qu = f( x，y) ，( x，y) ，xxyyBu = g( x，y) ，( x，y) ，其中 p 0，q 0 或者 p 0，q 0， R2 ，f( x，y) ，g( x，y) 是连续函数，B 是边界算子 在 内配置 N 个离散的

14、数据点 r1 ，r2 ，rN ，其中 r1 ，r2 ，rL 是内部节点，rL +1 ，rL +2 ，rN 是外部节点 设方程 ( 1)N的数值解为uN ( r)= j ( r rj ) ，( 5)j = 1其中 1 ，2 ，N 为待定系数 由 ( 1) 式和 ( 5) 式得Nj ( p( ri rj ) + q( ri rj ) )= f( ri ) ，i = 1，2，L( 6)j = 1Nj B( ri rj )= g( ri ) ，i = L + 1，L + 2，N，( 7)j = 1由 ( 6) 式，( 7) 式得到矩阵方程A = F， 其中( 8)( r1 r1 )( rL r1 )(

15、 r1 r2 )( rL r2 )B( rL +1 r2 )B( rN r2 )( r1 rN )( rL rN )æ ç ç çö÷÷÷A =ç÷ç B( rL +1 r1 )B( rL +1 rN ) ÷ç÷÷øB( rN rN )çè B( rN x1 )T ，FT =( 1 ，2 ，N )=( f( r1 ) ，f( rL ) ，g( rL +1 ) ，g( rN ) )由定理 2 知，矩阵方程 ( 8) 的

16、系数矩阵是非奇异的，所以矩阵方程有唯一解，只要从( 8) 式中解出 ，就N可以得到方程( 1) 的近似解 uN ( r) = j ( r rj ) j = 12数值算例在这一部分，我们将通过一个数值算例来讨论正定径向基函数解常微分方程，用绝对误差描述数值解与精确解的差别 绝对误差的形式如下，Error = maxE( x) Num( x)其中 E( x) 是精确解，Num( x)例 已知偏微分方程是数值解 ( uxx + uyy ) + 4u= 6sin xcos y， ( x，y) ，( x，y) ，u | = 0其中 = 0， × ， ， 是 的边界 此微分方程的精确解是 u=

17、sin xcos y22下面利用径向基函数配置法对上述微分方程求数值解，把区域 等分，步长 h ，分别选用 ( r)=42ce cr2 ( c 0，为形状参数) 和 ( r)差结果如表 1 和表 2( c 0，为形状参数) 作为基函数，所得数值解与精确解的绝对误=c2 + r2( r) = e cr2 为基函数，不同形状参数时的绝对误差的比较表 1hE( x)Error c = 0 03Error c = 0 05Error c = 0 5( ， )0 500 000 000 000 000 389 324 733 E 30 535 912 085 E 30 009 444 635 043 0

18、544( ，0)0 707 106 781 186 550 112 667 033 E 30 231 889 585 E 30 009 336 350 464 104 ， )(0 500 000 000 000 000 389 324 995 E 30 535 912 079 E 30 009 444 635 043 0544( ， )0 707 106 781 186 550 112 668 139 E 30 231 889 633 E 30 009 336 350 464 1024( ，0)1 000 000 000 000 000 103 084 777 E 30 136 895 125

19、 E 30 000 999 433 092 342 ， )(0 707 106 781 186 550 112 668 547 E 30 231 889 582 E 30 009 336 350 464 1024( 3， )0 500 000 000 000 000 389 325 141 E 30 535 912 032 E 30 009 444 635 043 0544( 3，0)0 707 106 781 186 550 112 668 430 E 30 231 889 545 E 30 009 336 350 464 104( 3， )0 500 000 000 000 000 389

20、 327 702 E 30 535 912 107 E 30 009 444 635 043 0544 )( 0，00 000 000 232 E 30 000 000 034 E 30 000 000 000 000 002( 0，0)00 000 000 261 E 30 000 000 157 E 30 000 000 000 000 00( 0， )00 000 002 095 E 30 000 000 145 E 30 000 000 000 000 002( ， )00 000 001 105 E 30 000 000 020 E 30 000 000 000 000 0022(

21、， )00 000 000 058 E 30 000 000 189 E 30 000 000 000 000 0022 )( ，00 000 001 629 E 30 000 000 125 E 30 000 000 000 000 002( ，0)00 000 002 270 E 30 000 000 114 E 30 000 000 000 000 00( ， )00 000 000 116 E 30 000 000 03 E 30 000 000 000 000 00 2 2c表 2 ( r)为基函数，不同形状参数时的绝对误差=c2 + r2hE( x)Error c = 7 6Err

22、or c = 0 05Error c = 0 5( ， )0 500 000 0000 011 132 791 042 430 496 193 059 036 290 006 269 771 246 3544( ，0)0 707 106 7810 017 381 333 564 240 701 736 878 297 950 025 874 233 158 374 ， )(0 500 000 0000 011 132 791 031 740 496 193 059 036 290 006 269 771 246 3544( ， )0 707 106 7810 017 381 333 569 4

23、70 701 736 878 297 950 025 874 233 158 3724( ，0)1 000 000 0000 025 526 035 609 570 992 421 578 520 720 028 027 307 902 722 ， )(0 707 106 7810 011 132 791 042 430 701 736 878 297 950 025 874 233 158 3724( 3， )0 500 000 0000 017 381 333 564 240 496 193 059 036 290 006 269 771 246 3544( 3，0)0 707 106 7

24、810 011 132 791 031 740 701 736 878 297 950 025 874 233 158 374( 3， )0 500 000 0000 017 381 333 569 470 496 193 059 036 290 006 269 771 246 3544 )( 0，00 025 526 035 609 570 000 000 000 000 000 000 000 000 000 002( 0，0)00 000 000 000 004 550 000 000 000 000 000 000 000 000 000 00( 0， )00 000 000 000

25、016 140 000 000 000 000 000 000 000 000 000 002( ， )00 000 000 000 005 230 000 000 000 000 000 000 000 000 000 0022( ， )00 000 000 000 002 730 000 000 000 000 000 000 000 000 000 0022 )( ，00 000 000 000 008 870 000 000 000 000 000 000 000 000 000 002( ，0)00 000 000 000 007 500 000 000 000 000 000 00

26、0 000 000 000 00( ， )00 000 000 000 010 230 000 000 000 000 000 000 000 000 000 00 2 表 1，表 2 取相同的步长 h ，但是径向基函数不同 表 1 以 ( r) = e cr2 为径向基函数，当取形状参=4数 c = 0 03 和 c = 0 05，绝对误差均达到 10 4 ，取形状参数 c = 0 5，最大绝对误差为 10 3 表 2 以径向基函2c 0，为形状参数) 作为基函数，当形状参数 c = 0 5，绝对误差达到 10 3 数 ( r)( c=c2 + r2通过表 1 和表 2 注意到，当取相同的形

27、状参数 c = 0 5，所得到的绝对误差都达到了 10 3 ，但是 c= 0 05时，表 2 的绝对误差明显不好 通过数值实验，当径向基函数 ( r) = e cr2 的形状参数 c 0 1 时，最大绝对2c误差为 10 3 ，当形状参数 c 0 1，最大绝对误差精确程度很差; 当径向基函数 ( r)=的形状参数c2 + r21 c 4 时，最大绝对误差为 0 1，当形状参数 4 c 10，最大绝对误差精确程度可以达到 10 3 ，而在 0 到1 之间取值时，变化比较大3结论本文成功地用正定径向基函数解一类偏微分方程，通过与精确解作比较，有很小的绝对误差，所以在插值点得到令人满意的数值解，并验

28、证了逼近过程中所得到的矩阵方程有唯一解在数值算例部分，选用了不同的径向基函数作为基函数进行数值逼近，当取合适的形状参数时，都可以 得到令人满意的数值解 但是注意到，选择相同的步长及不同的形状参数，有不同的绝对误差，所以所得数值 解与径向基函数形状参数的选取密切相关 如何选取形状参数，已经引起很多专家学者的注意 一般地，通过2 个形状参数，对所得数值解进行比较，逐步选取合适的形状参数，根据经验逐步计算，这种方法并不是很 好，甚至有时候不能得到满意的数值解 所以径向基函数的形状参数如何选取，还是需要进一步研究的问题另外，本文的数值实验部分，通过数值比较，用 Gauss 函数 ( r) = e cr

29、2 做为基函数明显优于用 ( r)=c2c2 + r2 作为基函数，那么，一类微分方程哪一类径向基函数可以得到比较好的数值解，这也是需要进一步研究的问题 一般来说，比较常用的是 Gauss 函数 ( r) = e cr2 作为基函数参考文献:123胡建伟，汤怀民 微分方程数值方法M 天津: 科学出版社，1999张池平 数值方法M 北京: 科学出版社，2006MADYCH W R Miscellaneous error bounds for multiquadric and related interpolatorsJ Comput Math Applic，1992，24( 12) : 121-

30、138KANSA E J Multiquadrics: a scarrered data approximation scheme with application to computational fluid dynamics: surface approximations and partial derivative estimatesJ Comput Math Applic，1990，19( 8 /9) : 127-154KANSA E J Multiquadrics: a scarrered data approximation scheme with application to c

31、omputational fluid dynamics: parabolic，and elliptic partial differential equationsJ Comput Math Applic，1990，19( 8 /9) : 146-161MADYCH W R，NELSON S A Miscellaneous error bounds for multiquadic and related inter polatorsJ Comput Math Appl，1992，24( 12) : 121-138MADYCH W R，NELSON S A，NELSON Multivariate interpolation and conditionally positive definite functionsJ Mat