matlab rbf函数_基于径向基函数（RBF）的无网格伪谱法与程序实现（2）—

参考资料

Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB.

P.387 P401

数值实现

Matlab 2019a

地球物理局 地震波动力学实验室 无网格组声明：
# 系列文章优先满足个人研究需求
# 欢迎批评指正，禁止转载

微分矩阵

为了了解如何找到一个微分矩阵，考虑展开式（1），并令

$B_{j}, j=1, \ldots, N$ 为一组任意线性无关的光滑函数，将作为我们近似空间的基。

如果我们在配点

$x_{i}, i=1, \ldots, N$ 处计算（1），则我们有

$\hat{u}\left(x_{i}\right)=\sum_{j=1}^{N} c_{j} B_{j}\left(x_{i}\right), \quad i=1, \ldots, N$

或者用矩阵-向量表示：

$\boldsymbol{u}={A} \mathbf{c}\quad\quad\quad(5)$

其中

$\boldsymbol{c}=\left[c_{1}, \ldots, {c}_{N}\right]^{T}$ 为系数向量，

评估矩阵（evaluation matrix）

$A$ 有项

$A_{i j}=B_{j}\left(x_{i}\right)$ ，

$\boldsymbol{u}$ 如前所述。

通过线性，我们还可以用展开式（1），通过对基函数进行微分，来计算

$\hat{{u}}$ ：

$\frac{d}{d x} \hat{u}(x)=\sum_{j=1}^{N} c_{j} \frac{d}{d x} B_{j}(x)$

如果我们再一次在配点

$x_{i}$ 评估，那么我们得到矩阵-向量表示：

$\boldsymbol{u}^{\prime}=A_{x} \boldsymbol{c}\quad\quad\quad(6)$

其中

$\boldsymbol{u}$ 和

$\boldsymbol{c}$ 与（5）中相同，而

导数矩阵（derivative matrix）

$A_{x}$ 有项

$\frac{d}{d x} B_{j}\left(x_{i}\right)$ ，或者，在径向基函数的情况下，

$\left.\frac{d}{d x} \varphi\left(\left\|x-x_{j}\right\|\right)\right|_{x=x_{i}}$ 。

为了获得微分矩阵

$D$ ，我们需要确保评估矩阵

$A$ 的可逆性。

这既取决于选择的基函数，也取决于配点

$x_{i}$ 的位置。对于单变量多项式，众所周知，对于任何不同的配点集，评估矩阵都是可逆的。特别是，如果多项式以基本（或拉格朗日）形式写出，则评估矩阵为单位矩阵。如果我们使用严格正定的径向基函数，则矩阵

$A$ 对于任何不同的配点集（也是非均匀间隔的点，且在

$\mathbb{R}^{s}, s>1$ 中）都是可逆的。另一方面，基本的径向基函数很难获得。对于统一的一维网格的特殊情况，可以在Platte and Driscoll(2005)中找到显式。在第14章中，对RBF的基本表示进行了一般性讨论。在下文中，我们将不坚持使用基本表示。

现在我们已经讨论了

$A$ 的可逆性，我们可以用（5）来正式求解系数向量

$\boldsymbol{c}=A^{-1} \boldsymbol{u}$ ，结合（6）得到：

$\boldsymbol{u}^{\prime}=A_{x} A^{-1} \boldsymbol{u}$

所以对应于（2）的微分矩阵

$D$ 由下式给出：

$D=A_{x} A^{-1}$

对于具有常系数的更复杂的线性微分算子

$\mathcal{L}$ ，我们以类似的方式进行操作，以获得离散化的微分算子（微分矩阵）：

$L=A_{\mathcal{L}} A^{-1}\quad\quad\quad(7)$

其中矩阵

$A_{\mathcal{L}}$ 有项

$\left(A_{\mathcal{L}}\right)_{i j}=\mathcal{L} B_{j}\left(x_{i}\right)$ 。在径向基函数的情况下，这些项具有形式

$\left(A_{\mathcal{L}}\right)_{i j}=\left.\mathcal{L} \varphi\left(\left\|x-x_{j}\right\|\right)\right|_{x=x_{i}}$ 。

在伪谱法的背景中，微分矩阵

$D$ 或

$L$ 现在可用于求解各种PDE（时间相关和时间无关）。例如，对于许多时间步进算法（如本章开头给出的例子），有时只需要乘以

$L$ 。对于其它问题，需要对

$L$ 求逆。在标准伪谱情况下，已知Chebyshev微分矩阵具有

$N$ 重零特征值（参见Canuto et al.(1998),p.70），因此，它本身是不可逆的。但是，一旦考虑到边界条件，情况就会改变（参见，例如，Trefethen(2000), p.67）。

例子

为了更深入地了解径向基函数的特殊性质，我们通过忽略边界条件，假定求解形式为

$\mathcal{L} u=f$ 的不适定线性椭圆PDE。通过求解离散线性方程组，可以获得在配点

$x_{i}$ 处的近似解：

$L \boldsymbol{u}=\boldsymbol{f}$

其中

$\boldsymbol{f}$ 包含了在配点处

$f$ 的值，而

$L$ 如前所述。换句话说，配点处的解由下式给出（参见(7)）：

$\boldsymbol{u}=L^{-1} \boldsymbol{f}=A\left(A_{\mathcal{L}}\right)^{-1} \boldsymbol{f}$

我们看到，为了进行下一步，

$L$ （以及

$A_{\mathcal{L}}$ ）的可逆性需要满足。

如上所述，基于Chebyshev多项式的伪谱法的微分矩阵是奇异的。这是很自然的，因为仅根据其导数的值重建未知函数的问题是不适定的。

但是，如果使用径向基函数，则引用广义Hermite插值结果可确保矩阵

$A_{\mathcal{L}}$ 是可逆的，前提是使用严格正定的基函数且微分算子为椭圆型。因此，基于RBF的伪谱法的基本微分矩阵

$L$ 是可逆的。

上面进行的观察表明，RBF方法有时“效果太好了以至于显得太假”。它们甚至可以为不适定的问题提供“解”。这是最优性原理的结果，即，由于本质空间范数的最小化变量，RBF方法具有内置的正则化功能。RBF的这一有趣功能最近已用于求解不适定问题（例如，参见Cheng and Cabral（2005））。

用MATLAB计算RBF微分矩阵

我们首先说明如何计算离散微分算子（微分矩阵）。

例如，为了计算一阶微分矩阵，我们需要记住——通过链式法则——RBF的导数将具有一般形式：

$\frac{\partial}{\partial x} \varphi(\|\boldsymbol{x}\|)=\frac{x}{r} \frac{d}{d r} \varphi(r)$

因此，我们不仅需要距离

$r$ ，还需要

$x$ 的微分，其中

$x$ 是

$\boldsymbol{x}$ 的第一个分量。因此，第一个MATLAB子函数中的主要语句是对这些距离和微分矩阵的计算。

    r = DistanceMatrix(x,x);        %   距离矩阵计算dx = DifferenceMatrix(x,x);     %   微分矩阵计算

DistanceMatrix.m

% Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB.
% 地球物理局 地震波动力学实验室 无网格组function DM = DistanceMatrix(dsites,ctrs)
%  输入：
%  dsites：M*s矩阵，代表一组R^s空间（即：每行包含一个s维点）中M数据点集
%  ctrs：N*s矩阵，代表一组R^s空间中N中心（每列一个中心）
%  DM：M*N矩阵，i,j位置包含第i个数据点和第j个中心之间的欧几里得距离[M,s] = size(dsites);[N,s] = size(ctrs);DM = zeros(M,N);
%  累积坐标差的平方和
%  ndgrid命令产生两个M*N矩阵：
%  dr：包含N个相同的列（每一列包含M数据点的第d个坐标）
%  cc：包含M个相同的行（每一行包含N中心的第d个坐标）for d = 1:s[dr,cc] = ndgrid(dsites(:,d),ctrs(:,d));   %   ndgrid：n维空间中的矩形网格DM = DM + (dr - cc).^2;endDM = sqrt(DM);
end

DifferenceMatrix.m

% Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB.
% P.342
%   地球物理局 地震波动力学实验室 无网格组function DM = DifferenceMatrix(datacoord,centercoord)
% 构成R中两组点的微分矩阵
% （R^s中有一些固定点坐标），即
% DM(j,k) = datacoord_j - centercoord_k
% ndgrid命令产生两个M*N矩阵
% dr和cc[dr,cc] = ndgrid(datacoord(:),centercoord(:));DM = dr - cc;
end

根据前面的讨论，微分矩阵由

$D=A_{x} A^{-1}$ 给出。这由下面语句实现：

    A = rbf(ep,r);Ax = dxrbf(ep,r,dx);D = Ax/A;

注意，在Matlab中使用

$/$ 或mrdivide来求解

$D A=A_{x}$ 。

为了完成程序，还要包括留一法交叉验证以优化RBF形状参数的步骤。这样对于矩阵问题

$D=A_{x} A^{-1} \Longleftrightarrow A^{T} D^{T}=\left(A_{x}\right)^{T}$ ，就可以最优化

$\varepsilon$ 的选择。

    mine = 0.1;    maxe = 10;              %   形参数区间ep = fminbnd(@(ep) CostEpsilonDRBF(ep,r,dx,rbf,dxrbf),mine,maxe);fprintf('Using epsilon = %fn',ep);

CostEpsilonDRBF.m

% Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB.
% P.402
%   地球物理局 地震波动力学实验室 无网格组function ceps = CostEpsilonDRBF(ep,r,dx,rbf,dxrbf)
%   提供epsilon的代价函数，用于形状参数LOOCV最优化
%   ep：形参数的值
%   r,dx：距离与微分矩阵
%   rbf,dxrbf：rbf及其导数的定义N = size(r,2);A = rbf(ep,r);      %   A^T 因为A是对称的rhs = dxrbf(ep,r,dx)';      %   A_x^TinvA = pinv(A);EF = (invA*rhs)./repmat(diag(invA),1,N);   % repmat：复制和平铺矩阵ceps = norm(EF(:));
end

现在我们计算右端的对应于矩阵

$A_{x}$ 的转置的矩阵。因此，需要通过repmat命令复制矩阵。

$\varepsilon$ 的代价现在是矩阵EF的Frobenius范数。针对该误差的其它度量也是合适的。对于标准插值设置，Rippa比较了

$\ell_{1}$ 和

$\ell_{2}$ 范数的使用（参见Rippa(1999)），并得出结论，

$\ell_{1}$ 范数产生更精确的“最优”。对于RBF-PS问题，我们观察到

$\ell_{2}$ 范数的结果非常好。

因此，综上所述，得到一维情况下计算微分矩阵的程序DRBF.m。

DRBF.m

% Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB.
% P.402
%   地球物理局 地震波动力学实验室 无网格组function [D,x] = DRBF(N,rbf,dxrbf)
%  计算对于一维导数的微分矩阵D
%  使用Chebyshev点，对于最优化形状参数采用LOOCV
%  输入：
%  N：创建N+1个配点
%  rbf,dxrbf：函数句柄，表示rbf及其导数if N == 0D = 0;x = 1;returnendx = cos(pi*(0:N)/N)';   %   Chebyshev点mine = 0.1;maxe = 10;              %   形参数区间r = DistanceMatrix(x,x);        %   距离矩阵计算dx = DifferenceMatrix(x,x);     %   微分矩阵计算ep = fminbnd(@(ep) CostEpsilonDRBF(ep,r,dx,rbf,dxrbf),mine,maxe);fprintf('Using epsilon = %fn',ep);A = rbf(ep,r);Ax = dxrbf(ep,r,dx);D = Ax/A;
end