參考資料

Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB.

P.387 P401

數(shù)值實(shí)現(xiàn)

Matlab 2019a

地球物理局 地震波動(dòng)力學(xué)實(shí)驗(yàn)室 無(wú)網(wǎng)格組聲明： # 系列文章優(yōu)先滿足個(gè)人研究需求 # 歡迎批評(píng)指正，禁止轉(zhuǎn)載

目錄

石中居士：基于徑向基函數(shù)（RBF）的無(wú)網(wǎng)格偽譜法與程序?qū)崿F(xiàn)——目錄?zhuanlan.zhihu.com

---

微分矩陣

為了了解如何找到一個(gè)微分矩陣，考慮展開式（1），并令

為一組任意線性無(wú)關(guān)的光滑函數(shù)，將作為我們近似空間的基。

如果我們?cè)谂潼c(diǎn)

處計(jì)算（1），則我們有

或者用矩陣-向量表示：

其中

為系數(shù)向量，評(píng)估矩陣（evaluation matrix） 有項(xiàng) ， 如前所述。

通過(guò)線性，我們還可以用展開式（1），通過(guò)對(duì)基函數(shù)進(jìn)行微分，來(lái)計(jì)算

：

如果我們?cè)僖淮卧谂潼c(diǎn)

評(píng)估，那么我們得到矩陣-向量表示：

其中

和 與（5）中相同，而導(dǎo)數(shù)矩陣（derivative matrix） 有項(xiàng) ，或者，在徑向基函數(shù)的情況下， 。

為了獲得微分矩陣

，我們需要確保評(píng)估矩陣 的可逆性。這既取決于選擇的基函數(shù)，也取決于配點(diǎn) 的位置。對(duì)于單變量多項(xiàng)式，眾所周知，對(duì)于任何不同的配點(diǎn)集，評(píng)估矩陣都是可逆的。特別是，如果多項(xiàng)式以基本（或拉格朗日）形式寫出，則評(píng)估矩陣為單位矩陣。如果我們使用嚴(yán)格正定的徑向基函數(shù)，則矩陣 對(duì)于任何不同的配點(diǎn)集（也是非均勻間隔的點(diǎn)，且在 中）都是可逆的。另一方面，基本的徑向基函數(shù)很難獲得。對(duì)于統(tǒng)一的一維網(wǎng)格的特殊情況，可以在Platte and Driscoll(2005)中找到顯式。在第14章中，對(duì)RBF的基本表示進(jìn)行了一般性討論。在下文中，我們將不堅(jiān)持使用基本表示。

現(xiàn)在我們已經(jīng)討論了

的可逆性，我們可以用（5）來(lái)正式求解系數(shù)向量 ，結(jié)合（6）得到：

所以對(duì)應(yīng)于（2）的微分矩陣

由下式給出：

對(duì)于具有常系數(shù)的更復(fù)雜的線性微分算子

，我們以類似的方式進(jìn)行操作，以獲得離散化的微分算子（微分矩陣）：

其中矩陣

有項(xiàng) 。在徑向基函數(shù)的情況下，這些項(xiàng)具有形式 。

在偽譜法的背景中，微分矩陣

或 現(xiàn)在可用于求解各種PDE（時(shí)間相關(guān)和時(shí)間無(wú)關(guān)）。例如，對(duì)于許多時(shí)間步進(jìn)算法（如本章開頭給出的例子），有時(shí)只需要乘以 。對(duì)于其它問(wèn)題，需要對(duì) 求逆。在標(biāo)準(zhǔn)偽譜情況下，已知Chebyshev微分矩陣具有 重零特征值（參見(jiàn)Canuto et al.(1998),p.70），因此，它本身是不可逆的。但是，一旦考慮到邊界條件，情況就會(huì)改變（參見(jiàn)，例如，Trefethen(2000), p.67）。

例子

為了更深入地了解徑向基函數(shù)的特殊性質(zhì)，我們通過(guò)忽略邊界條件，假定求解形式為

的不適定線性橢圓PDE。通過(guò)求解離散線性方程組，可以獲得在配點(diǎn) 處的近似解：

其中

包含了在配點(diǎn)處 的值，而 如前所述。換句話說(shuō)，配點(diǎn)處的解由下式給出（參見(jiàn)(7)）：

我們看到，為了進(jìn)行下一步，

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

如上所述，基于Chebyshev多項(xiàng)式的偽譜法的微分矩陣是奇異的。這是很自然的，因?yàn)閮H根據(jù)其導(dǎo)數(shù)的值重建未知函數(shù)的問(wèn)題是不適定的。

但是，如果使用徑向基函數(shù)，則引用廣義Hermite插值結(jié)果可確保矩陣

是可逆的，前提是使用嚴(yán)格正定的基函數(shù)且微分算子為橢圓型。因此，基于RBF的偽譜法的基本微分矩陣 是可逆的。

上面進(jìn)行的觀察表明，RBF方法有時(shí)“效果太好了以至于顯得太假”。它們甚至可以為不適定的問(wèn)題提供“解”。這是最優(yōu)性原理的結(jié)果，即，由于本質(zhì)空間范數(shù)的最小化變量，RBF方法具有內(nèi)置的正則化功能。RBF的這一有趣功能最近已用于求解不適定問(wèn)題（例如，參見(jiàn)Cheng and Cabral（2005））。

---

用MATLAB計(jì)算RBF微分矩陣

我們首先說(shuō)明如何計(jì)算離散微分算子（微分矩陣）。

例如，為了計(jì)算一階微分矩陣，我們需要記住——通過(guò)鏈?zhǔn)椒▌t——RBF的導(dǎo)數(shù)將具有一般形式：

因此，我們不僅需要距離

，還需要 的微分，其中 是 的第一個(gè)分量。因此，第一個(gè)MATLAB子函數(shù)中的主要語(yǔ)句是對(duì)這些距離和微分矩陣的計(jì)算。 r = DistanceMatrix(x,x); % 距離矩陣計(jì)算dx = DifferenceMatrix(x,x); % 微分矩陣計(jì)算

DistanceMatrix.m

% Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB. % 地球物理局 地震波動(dòng)力學(xué)實(shí)驗(yàn)室 無(wú)網(wǎng)格組function DM = DistanceMatrix(dsites,ctrs) % 輸入： % dsites：M*s矩陣，代表一組R^s空間（即：每行包含一個(gè)s維點(diǎn)）中M數(shù)據(jù)點(diǎn)集 % ctrs：N*s矩陣，代表一組R^s空間中N中心（每列一個(gè)中心） % DM：M*N矩陣，i,j位置包含第i個(gè)數(shù)據(jù)點(diǎn)和第j個(gè)中心之間的歐幾里得距離[M,s] = size(dsites);[N,s] = size(ctrs);DM = zeros(M,N); % 累積坐標(biāo)差的平方和 % ndgrid命令產(chǎn)生兩個(gè)M*N矩陣： % dr：包含N個(gè)相同的列（每一列包含M數(shù)據(jù)點(diǎn)的第d個(gè)坐標(biāo)） % cc：包含M個(gè)相同的行（每一行包含N中心的第d個(gè)坐標(biāo)）for d = 1:s[dr,cc] = ndgrid(dsites(:,d),ctrs(:,d)); % ndgrid：n維空間中的矩形網(wǎng)格DM = DM + (dr - cc).^2;endDM = sqrt(DM); end

DifferenceMatrix.m

% Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB. % P.342 % 地球物理局 地震波動(dòng)力學(xué)實(shí)驗(yàn)室 無(wú)網(wǎng)格組function DM = DifferenceMatrix(datacoord,centercoord) % 構(gòu)成R中兩組點(diǎn)的微分矩陣 % （R^s中有一些固定點(diǎn)坐標(biāo)），即 % DM(j,k) = datacoord_j - centercoord_k % ndgrid命令產(chǎn)生兩個(gè)M*N矩陣 % dr和cc[dr,cc] = ndgrid(datacoord(:),centercoord(:));DM = dr - cc; end

根據(jù)前面的討論，微分矩陣由

給出。這由下面語(yǔ)句實(shí)現(xiàn)： A = rbf(ep,r);Ax = dxrbf(ep,r,dx);D = Ax/A;

注意，在Matlab中使用

或mrdivide來(lái)求解 。

為了完成程序，還要包括留一法交叉驗(yàn)證以優(yōu)化RBF形狀參數(shù)的步驟。這樣對(duì)于矩陣問(wèn)題

，就可以最優(yōu)化 的選擇。 mine = 0.1; maxe = 10; % 形參數(shù)區(qū)間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 % 地球物理局 地震波動(dòng)力學(xué)實(shí)驗(yàn)室 無(wú)網(wǎng)格組function ceps = CostEpsilonDRBF(ep,r,dx,rbf,dxrbf) % 提供epsilon的代價(jià)函數(shù)，用于形狀參數(shù)LOOCV最優(yōu)化 % ep：形參數(shù)的值 % r,dx：距離與微分矩陣 % rbf,dxrbf：rbf及其導(dǎo)數(shù)的定義N = size(r,2);A = rbf(ep,r); % A^T 因?yàn)锳是對(duì)稱的rhs = dxrbf(ep,r,dx)'; % A_x^TinvA = pinv(A);EF = (invA*rhs)./repmat(diag(invA),1,N); % repmat：復(fù)制和平鋪矩陣ceps = norm(EF(:)); end

現(xiàn)在我們計(jì)算右端的對(duì)應(yīng)于矩陣

的轉(zhuǎn)置的矩陣。因此，需要通過(guò)repmat命令復(fù)制矩陣。 的代價(jià)現(xiàn)在是矩陣EF的Frobenius范數(shù)。針對(duì)該誤差的其它度量也是合適的。對(duì)于標(biāo)準(zhǔn)插值設(shè)置，Rippa比較了 和 范數(shù)的使用（參見(jiàn)Rippa(1999)），并得出結(jié)論， 范數(shù)產(chǎn)生更精確的“最優(yōu)”。對(duì)于RBF-PS問(wèn)題，我們觀察到 范數(shù)的結(jié)果非常好。

因此，綜上所述，得到一維情況下計(jì)算微分矩陣的程序DRBF.m。

DRBF.m

% Gregory E. Fasshauer. Meshfree Approximation Methods with MATLAB. % P.402 % 地球物理局 地震波動(dòng)力學(xué)實(shí)驗(yàn)室 無(wú)網(wǎng)格組function [D,x] = DRBF(N,rbf,dxrbf) % 計(jì)算對(duì)于一維導(dǎo)數(shù)的微分矩陣D % 使用Chebyshev點(diǎn)，對(duì)于最優(yōu)化形狀參數(shù)采用LOOCV % 輸入： % N：創(chuàng)建N+1個(gè)配點(diǎn) % rbf,dxrbf：函數(shù)句柄，表示rbf及其導(dǎo)數(shù)if N == 0D = 0;x = 1;returnendx = cos(pi*(0:N)/N)'; % Chebyshev點(diǎn)mine = 0.1;maxe = 10; % 形參數(shù)區(qū)間r = DistanceMatrix(x,x); % 距離矩陣計(jì)算dx = DifferenceMatrix(x,x); % 微分矩陣計(jì)算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

---

總結(jié)

以上是生活随笔為你收集整理的matlab rbf函数_基于径向基函数（RBF）的无网格伪谱法与程序实现（2）——微分矩阵...的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。