0 引 言 Helmholtz方程作为简化的波动方程，长期以来备受关注.求解Helmholtz方程的数值方法主要有有限差分法、有限元法和边界元法等.这些方法主要基于网格近似或者需要背景网格积分，在处理高波数或高频率Helmholtz方程及高维问题时存在网格划分困难和内存消耗大等缺点.基于点的近似思想，无网格方法在一定程度上减少对网格的依赖，甚至不需要网格.目前，主要有2类无网格法：区域型无网格法和边界型无网格法.区域型无网格法主要有Kansa方法和MLPG法等；边界型无网格法由于只需要边界节点的信息，近年来备受关注.基本解法作为一个典型的边界型无网格法已被成功应用于工程领域，如数值求解Helmholtz方程等，基本解法已经取得很大的成功，但为克服核函数的奇异性，虚拟边界的选取至今仍无定论.CHEN等［1-3］提出的边界节点法不仅保持基本解法的优点，而且克服基本解法最大的缺点――虚拟边界难以选取.边界节点法的主要思想是将控制方程的通解作为核函数，避免基本解法中核函数的奇异性，将源点和配点取在物理边界上.边界节点法已被成功应用于Helmholtz�方程求解.
为使研究人员更好地了解和运用边界节点法，本文开发基于边界节点法的MATLAB工具箱.自MATLAB出现以来，各领域的专家为研究高技术和应用方面的问题，已经开发出符号计算工具箱和图形处理工具箱等各领域的工具箱，而边界节点法工具箱仍未见报道.本文介绍基于边界节点法开发的MATLAB工具箱，用于求解简化的Helmholtz方程，主要考虑齐次Helmholtz方程、非齐次Helmholtz方程以及Helmholtz方程控制的边界型反问题.运用本文开发的MATLAB工具箱，可非常简单地运用边界节点法拟合上述问题，将正则化方法加入到本文开发的MATLAB工具箱中，便于研究人员更好地了解正则化方法对边界节点法的影响.
1 Helmholtz方程和边界节点法
1.1 Helmholtz方程
二维Helmholtz控制方程及其边界条件为式中：
Δ为Laplace算子；
�λ�为控制方程的波数；
�Ω�为二维单连通区域；
�Γ�1和Γ��2分别为Dirichlet边界和Neumann边界；
�B�1和B��2为边界条件算子；
�p�为二维区域节点坐标�(x,y)�；
�f(p)�为方程的源项；
�g(p�1)和�h(p�2)�分别为给定的Dirichlet和Neumann边界条件值.当f(p)=0时，式（1）为齐次方程，反之为非齐次方程.当Γ�1∩Γ�2=�且Γ�1∪Γ�2=Ω时，只需要给出区域内部的数值解；当同时给出边界Γ�1上的�Dirichlet条件和Neumann�条件且边界Γ�2不可测时，式（1）～（3）构成一类典型的反问题，需给出不可测边界Γ�2上的数值解.
1.2 边界节点法
边界节点法主要通过控制方程的非奇异通解的线性组合来近似方程的解，可表示为当N增加到一定程度时，式（6）的条件数会剧烈变化，所得解的收敛曲线会出现振荡现象.有效条件数能更好地表现矩阵的病态程度，故引入有效条件数.当求解的问题为反问题时，矩阵方程的解并不连续依赖于边界条件的值，且若边界条件含有扰动时，会对数值解带来致命影响.正则化方法已经被用于解决上述问题.［4-5］本文介绍的MATLAB工具箱主要采用以下几种正则化方法.
2 正则化方法和正则化参数
在引入正则化方法前简单介绍奇异值分解（Singular Value Decomposition,SVD）.对于式（6）中的矩阵式中：第一项为残差项；第二项为正则化项，包含解的先验信息，如光滑性和有界性等；参数�μ�≥0为正则化参数控制正则化项和残差项的比重，当�μ�=0时，式(16)转化为最小二乘问题.基于SVD的Tikhonov正则化方法可描述为
3 MATLAB工具箱简介
基于边界节点法求解Helmholtz方程的MATLAB工具箱主要由3类功能型的函数组成：(1)产生边界数据及内部检测点数据函数;(2)公用数据计算函数;(3)求解函数.
3.1 产生边界数据和内部检测点数据函数
3.1.1 产生边界数据函数
bkm_2d_circle
功能：产生二维圆形边界数据.
用法：�［x，y］�=bkm_2d_circle(�O�x，O�y，r，n�).
描述：用于产生以�（O�x，O�y）�为圆心、�r�为半径的均匀分布的�n�个点.
bkm _2d_rect
功能：用于产生二维矩形边界节点.
用法：�［x，y］�=bkm_2d_rect(�O�x，O�y，a，b，an，bn�).
描述：用于在（�O�x，O�y�）为中心、长为�a宽为b�的矩形的边界上产生边界点，其中两边上
均匀分布着�an和bn�个点.
bkm _2d_ellip
功能：用于产生二维椭圆形的边界数据.
用法：［�x，y�］= bkm _2d_ellip(�O�x，O�y，a，b，n�).
描述：用于产生椭圆形的数据，其中椭圆的中心为（�O�x，O�y�），长短轴分别为�a和b�
，边界节点数为�n�.
3.1.2 产生内部检测点函数
bkm_2d_circle_in
功能：用于产生圆内检测点.
用法：［�x，y�］= bkm_2d_circle_in(�O�x，O�y，r，n�，“method”).
描述：用于产生位于以（�O�x，O�y�）为圆心，�r�为半径的圆内分布的�n�个检测点.m
ethod=regular表示有规律分布的检测点；method=irregular表示无规律分布的检测点.
bkm _2d_rect_in
功能：用于产生矩形内部检测点.
用法：［�x，y�］= bkm _2d_rect_in（�O�x，O�y，a，b，n�，“method”）.
描述：用于产生以（�O�x，O�y�）为中心，边长为�a和b�的矩形内部检测点，个数为�n
�.method=regular表示有规律分布的检测点；method=irregular表示无规律分布的检测
点.
bkm _2d_ellip _in
功能：用于产生椭圆形内部检测点.
用法：［�x，y�］= bkm _2d_ellip _in（�O�x，O�y，a，b，n�，“method”）.
描述：用于产生以（�O�x，O�y�）为中心，长短轴长分别为�a和b�的椭圆内部分布的�n
�个点.method=regular表示有规律分布的检测点；method=irregular表示无规律分布的
检测点.
3.2 公用数据计算函数
bkm _2d_ rmatrix
功能：用于计算2组节点数据之间的距离矩阵.
用法：［�rM，vx，vy�］= bkm _2d_ rmatrix�(x�1，y�1，x�2，y�2�).
描述：用于计算数据（x�1，y�1）与数据（x�2，y�2）之间的距离矩阵.
rM(i，j)=(x�1(i)-x�2(j))�2+(y�1(i)-y�2(j))�2.其中，i为x�1
的长度，j为x�2的长度.
bkm _2d _cancl
功能：用于计算表达式的值，可用于求精确解或边界条件.
用法：［result］= bkm _2d _cancl(�x，y�，“method1”,“method2”，tip).
描述：用于计算表达式method1和method2的值，其中，�x，y�及tip为同维数的矩阵，当
tip(�i，j�)=1时计算表达式method1的值，否则计算表达式method2的值.
result(�i，j�)=method(�x(i，j)，y(i，j�)).
solveab
功能：用于计算矩阵方程
ax=
b，求解
x.
用法：�alpha=solveab�(
a，
b，“method”).
描述：用于计算矩阵方程�ax=b，其中，a�为系数矩阵，b[WT]为方程右端项.method=�A＼b�则直接调用MATLAB内部函数.若method=regu则要求输入正则化方法和正则化参数，可供选择的正则化方法有“TR”，“TSVD”和“DSVD”，可供选择的正则化参数有“L”和“GCV”.通过这些方法，最终可求出方程的未知数alpha.effCond
功能：用于计算矩阵的有效条件数.
用法：result=effCond(�a，b�).
描述：用于计算矩阵方程�ax=b�的有效条件数.
cond
功能：MATLAB自带用于求解矩阵条件数.
用法：result=cond(�a�).
描述：用于计算矩阵�a�的条件数.
3.3 求解函数
bkm_2d_helm
功能：用于求解二维Helmholtz齐次正问题和反问题，可得到矩阵方程的条件数、有效条件数、数值解、相对误差和绝对误差等.
用法：［cond，econd，appro，AerL，RerL］=bkm_2d_helm(lamda，�x，y�，result_�xy�，tip1，“method”， �x���in，�y���in， tip2， result_�xy���in).
描述：用于求解齐次正问题或反问题.其中，lamda为控制方程的波数；�x和y�分别为边界节点坐标；tip1为边界条件类型.当tip1（�i，j�）=1时，节点（�x(i，j)，y(i，j�)）为Dirichlet条件，否则为Neumann边界条件.result_�xy�为边界条件值.method=�A�＼b[WT]”或者method=regu用于标示求解矩阵方程的方法.当求解问题为齐次正问题时，
�x���in和�y���in为内部检测节点的坐标；当求解问题为反问题时，�x���in和�y���in为求解问题的不可测边界节点坐标.result_�xy���in为�x���in和�y���in上用于比较的精确解.tip2标示�x���in和�y���in的边界条件类型.cond及econd分别为求解的矩阵方程条件数和有效条件数，appro为(�x���in，�y���in)上的数值解.AerL为绝对误差，RerL为相对误差.
bkm_2d_helm_inho功能：用于求解二维非齐次Helmholtz方程.
用法：［cond，econd，appro，AerL，RerL］=bkm_2d_helm_inho (lamda，�x，y�，result_�xy�，tip1，“method”， �x���in，�y���in，tip2，result_�xy���in，�f�).
描述：用于求解非齐次正问题或反问题.其中，lamda为控制方程的波数；�x和y�分别为边界节点坐标；tip1为边界条件类型.当tip1（�i，j�）=1时，节点（�x(i，j)，y(i，j�)）为Dirichlet条件，否则为Neumann边界条件.result_�xy�为边界条件值.method=“�A�＼b[WT]”或者method=regu用于标示求解矩阵方程的方法.当求解问题为非齐次正问题时，�x���in和�y���in为内部检测节点的坐标.当求解问题为反问题时，�x���in和�y���in为求解问题的不可测边界节点坐标.result_�xy���in为�x���in和�y���in上用于比较的精确解.tip2标示�x���in和�y���in的边界条件类型.�f�为边界节点的源项，与�x和y�同维数.cond及econd分别为求解矩阵方程的条件数和有效条件数，appro为(�x���in，�y���in)上的数值解.AerL为绝对误差，RerL为相对误差.
4 算 例
算例1 本算例考虑齐次Helmholtz方程.求解区域为以（0，0）为中心、半径为1的圆，用于比较的解析解为�u（x，y）�=sin �x� cos �y�.配点为均匀分布在该圆上的10个点，其中，前8个点给出Dirichlet边界条件,后2个点给出Neumann边界条件.检测点规律分布在该圆内，总数为20.MATLAB求解过程为
>>［�x,y�］= bkm_2d_circle (0,0,1,10);
>>result_�xy�= bkm _2d _cancl (�x,y�,‘sin �x� cos �y�’，
‘(cos �x� cos �y�)�x�-(sin �x� sin �y�)�y�’,
［1 1 1 1 1 1 1 1 0 0］);
>>［�x���in,�y���in］=bkm_2d_circle _in(0,0,1,20,‘regular
’);
>>result_�xy���in=bkm_2d_cancl(�x���in,�y���in,�‘sin �x�cos �y�’，‘(cos �x� cos �y�)�x�-(sin �x� sin �y�)
�y�’，ones(size(�x���in)));
>>［cond,econd,appro_�xy���in,AerL,RerL］=
�bkm_2d_helm(sqrt(2),�x,y�，result_�xy�, ［1 1 1 1 1 1 1 1 0 0］， ‘�A�＼
b[WTBZ]’，
�x���in，�y���in，ones(size(�x���in,1)，size(�x���in，2)
),result_�xy���in)
得到的结果为
cond =8.203 3E+004
econd=218.153 0
appro_�xy���in=
0.098 4 -0.098 3 -0.198 7
-0.098 3 0.098 4 0.186 9
-0.186 9 -0.389 4 -0.186 8
0.186 9 0.256 5 -0.256 5
-0.564 6 -0.256 5 0.256 6
0.299 7 -0.299 6 -0.717 3
-0.299 6 0.299 7
AerL =2.817 9E-005
RerL =1.258 3E-004
为说明此算法的有效性，计算出使用基本解方法［5-6］所得到的结果.其中，边界节点值与算例1中边界节点法取值相同；虚拟边界选取在与此求解区域同心的圆上，半径为1.5；辅助点为均匀分布在辅助边界上的10个点.计算结果为AerL=0.076 5，RerL=0.379 9，可知，计算精度远小于边界节点法的精度，且边界节点法不受虚拟边界的约束.算例2 本算例考虑非齐次Helmholtz方程.求解区域为以（0,0）为中心、半径为1的圆，用于比较的解析解为�u（x，y�）=sin �x� sin �y�，非齐次项��f(x，y)�=0.配点为均匀分布在该圆上的10个点，其中，前8个点给出Dirichlet边界条件、后2个点给出Neumann边界条件.检测点规律分布在该圆内,总数为20.
MATLAB求解过程为
>>［�x，y�］= bkm_2d_circle (0,0,1,10);
>>result_�xy�= bkm _2d _cancl (�x,y�,‘sin �x� sin �y�)’，
‘(cos �x� sin �y�)・�x�+(sin �x� cos �y�)・�y�’，
［1 1 1 1 1 1 1 1 0 0］);
>>［�x���in,�y���in］= bkm_2d_circle _in(0,0,1,20，‘regula
r’);
>>result_�xy���in=bkm_2d_cancl(�x���in,�y���in，
‘sin �x�sin �y�’，‘(cos �x� sin �y�)・�x�+(sin �x� cos �y�)・
�y�’，ones(size(�x���in)));
>>［cond,econd,appro_�xy���in,AerL,RerL］=
�bkm_2d_helm(sqrt(2)，�x,y�,result_�xy�， ［1 1 1 1 1 1 1 1 0 0］，
‘regu’，�x���in，�y���in，
ones(size(�x���in,1),size(�x���in,2)),result_�xy���in)
得到的结果为
Please input the regularization method: ‘DSVD’ or ‘TSVD’ or ‘TR’ and the
regularization parameter: ‘L’ or ‘GCV’
Please input the regularization method:
‘DSVD’
Please input the regularization parameter:
‘L’
cond=8.203 3E+004
econd=7.039 5E+003
appro_�xy���in=-0.021 5 -0.005 4 -0.072 2
0.327 3 -0.112 9 -0.282 5
-0.212 8 …
AerL =2.008 7E-004
RerL =0.003 1
为说明算例2的有效性，计算出使用基本解方法所得的结果.其中，边界节点值与此算例中边界节点法中取值相同；虚拟边界选取在与此求解区域同心的圆上，半径为1.5；辅助点为均匀分布在辅助边界上的10个点.计算结果为AerL=0.035 7，RerL=0.543 5，同样可知，由基本解法所得到的结果精度远小于边界节点法的精度，且边界节点法不受虚拟边界的约束.
5 结束语
基于边界节点法的MATLAB工具箱具有使用方便、计算精度高等优点，通过选取适当的MATLAB函数，可将边界节点法用于求解声学问题.
参考文献：
�［1］ CHEN Wen, TANAKA M. A meshless, integration-free, and boundary-only RBF technique［J］. Computers & Math Applications, 2002, �43(3-5): 379-391.
［2］ CHEN Wen. Symmetric boundary knot method［J］. Eng Anal Boundary Elements, 2002, 26(6): 489-494.
［3］ CHEN Wen, TANAKA M. New insights into boundary-only and domain-type RBF methods［J］. Int J Nonlinear Sci & Numer Simulation, 2000, �1(3): 145-151.
［4］ HANSEN P C. Regularization tools: a MATLAB package for analysis and solution of discrete ill-posed problems［J］. Numer Algorithms, 1994, �6(1): 1-35.
［5］ WEI Ting, HON Y C, LING Leevan. Method of fundamental solutions with regularization techniques for Cauchy problems of elliptic operators［J］. Eng Anal Boundary Elements, 2007, 31(4): 373-385.
［6］ FAIRWEATHER G, KARAGEORGHIS A. The method of fundamental solutions for elliptic boundary value problems［J］. Adv Comput Math, 1998, 9(1): 69-95.
(编辑 陈锋杰)