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); endDifferenceMatrix.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）——微分矩阵...的全部內容，希望文章能夠幫你解決所遇到的問題。
