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

参考资料

Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB.

P.387 P401

数值实现

Matlab 2019a

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

目录

石中居士：基于径向基函数（RBF）的无网格伪谱法与程序实现——目录​zhuanlan.zhihu.com

---

微分矩阵

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

为一组任意线性无关的光滑函数，将作为我们近似空间的基。

如果我们在配点

处计算（1），则我们有

或者用矩阵-向量表示：

其中

为系数向量，评估矩阵（evaluation matrix） 有项 ， 如前所述。

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

：

如果我们再一次在配点

评估，那么我们得到矩阵-向量表示：

其中

和 与（5）中相同，而导数矩阵（derivative matrix） 有项 ，或者，在径向基函数的情况下， 。

为了获得微分矩阵

，我们需要确保评估矩阵 的可逆性。这既取决于选择的基函数，也取决于配点 的位置。对于单变量多项式，众所周知，对于任何不同的配点集，评估矩阵都是可逆的。特别是，如果多项式以基本（或拉格朗日）形式写出，则评估矩阵为单位矩阵。如果我们使用严格正定的径向基函数，则矩阵 对于任何不同的配点集（也是非均匀间隔的点，且在 中）都是可逆的。另一方面，基本的径向基函数很难获得。对于统一的一维网格的特殊情况，可以在Platte and Driscoll(2005)中找到显式。在第14章中，对RBF的基本表示进行了一般性讨论。在下文中，我们将不坚持使用基本表示。

现在我们已经讨论了

的可逆性，我们可以用（5）来正式求解系数向量 ，结合（6）得到：

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

由下式给出：

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

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

其中矩阵

有项 。在径向基函数的情况下，这些项具有形式 。

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

或 现在可用于求解各种PDE（时间相关和时间无关）。例如，对于许多时间步进算法（如本章开头给出的例子），有时只需要乘以 。对于其它问题，需要对 求逆。在标准伪谱情况下，已知Chebyshev微分矩阵具有 重零特征值（参见Canuto et al.(1998),p.70），因此，它本身是不可逆的。但是，一旦考虑到边界条件，情况就会改变（参见，例如，Trefethen(2000), p.67）。

例子

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

的不适定线性椭圆PDE。通过求解离散线性方程组，可以获得在配点 处的近似解：

其中

包含了在配点处 的值，而 如前所述。换句话说，配点处的解由下式给出（参见(7)）：

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

（以及 ）的可逆性需要满足。

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

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

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

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

---

用MATLAB计算RBF微分矩阵

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

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

因此，我们不仅需要距离

，还需要 的微分，其中 是 的第一个分量。因此，第一个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

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

给出。这由下面语句实现： A = rbf(ep,r);Ax = dxrbf(ep,r,dx);D = Ax/A;

注意，在Matlab中使用

或mrdivide来求解 。

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

，就可以最优化 的选择。 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

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

的转置的矩阵。因此，需要通过repmat命令复制矩阵。 的代价现在是矩阵EF的Frobenius范数。针对该误差的其它度量也是合适的。对于标准插值设置，Rippa比较了 和 范数的使用（参见Rippa(1999)），并得出结论， 范数产生更精确的“最优”。对于RBF-PS问题，我们观察到 范数的结果非常好。

因此，综上所述，得到一维情况下计算微分矩阵的程序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

---

总结

以上是生活随笔为你收集整理的matlab rbf函数_基于径向基函数（RBF）的无网格伪谱法与程序实现（2）——微分矩阵...的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。