功能:采用径向紧支基函数进行二维插值(Radial Basis Function)

格式:
InterpolationRBF2D(x0,y0,f0,x,y,n,a,k,t)
InterpolationRBF2D(x0,y0,f0,x,y,n,a,k)
InterpolationRBF2D(x0,y0,f0,x,y,n,a)
InterpolationRBF2D(x0,y0,f0,x,y,n)
InterpolationRBF2D(x0,y0,f0,x,y)

x0:矩阵变量存储的离散数据x坐标点
y0:矩阵变量存储的离散数据y坐标点,此参数与x0变量矩阵维数一致
f0:矩阵变量存储的离散数据在(x,y)坐标下对应的函数值,,此参数与x0、y0变量矩阵维数一致
x:矩阵变量存储的需要插值点的x坐标
y:矩阵变量存储的需要插值点的y坐标,此参数与x变量矩阵维数一致
n:选取离插值点最近的已知点个数,默认为4
a:半径动态控制倍数参数(>=1),默认为16
k:一个1*2的矩阵变量,其元素分别表示x0与y0对插值f0的权重系数,默认均一样,均取1。注意,权重系数不能小于0.
t:紧支基函数类型(此值从1到10)，具体对应函数可参看原理。默认为5。


说明:
函数成功返回与x、y变量同维数的插值矩阵。如果插值结果不理想,可尝试更改相应参数。最后的数据如果平滑度不够,可以使用函数{SmoothGrid<矩阵运算\SmoothGrid>}对结果进行加工处理。

原理:

1、对x点进行插值时,找出与x临近的n个点x0,采用本文算法计算各个半径比r

2、采用径向基函数进行求解

3、径向基函数类型如下(t对应的值)
(1)t=1-10来自文献[1]
t=1,  f(r)=1-r
t=2,  f(r)=(1-r)^3*(3*r+1)
t=3,  f(r)=(1-r)^5*(8*r^2+5*r+1)
t=4,  f(r)=(1-r)^2
t=5,  f(r)=(1-r)^4*(4*r+1)
t=6,  f(r)=(1-r)^6*(35*r^2+18*r+3)
t=7,  f(r)=(1-r)^8*(32*r^3+25*r^2+8*r+1)
t=8,  f(r)=(1-r)^3
t=9,  f(r)=(1-r)^5*(5*r+1)
t=10, f(r)=(1-r)^7*(16*r^2+7*r+1)

(2)t=11-19来自文献[2]
t=11, f(r)=(1-r)^2*(2+r)
t=12, f(r)=(1-r)^3*(1+3*r+r^2)
t=13, f(r)=(1-r)^3*(8+9*r+3*r^2)
t=14, f(r)=(1-r)^4*(4+16*r+12*r^2+3*r^3)
t=15, f(r)=(1-r)^5*(1+5*r+9*r^2+5*r^3+r^4)
t=16, f(r)=(1-r)^4*(16+29*r+20*r^2+5*r^3)
t=17, f(r)=(1-r)^5*(8+40*r+48*r^2+25*r^3+5*r^4)
t=18, f(r)=(1-r)^6*(6+36*r+82*r^2+72*r^3+30*r^4+5*r^5)
t=19, f(r)=(1-r)^7*(5+35*r+101*r^2+147*r^3+101*r^4+35*r^5+5*r^6)

参考:
/*
[1] Wendland H.Piecewise polynomial,positive definite and compactly supported radial functions of minimal degree[J]. Advances in computational Mathematics,1995,4(1):389-396. 

[2] Wu Z. Compactly supported positive definite radial functions[J]. Advances in Computational Mathematics, 1995, 4(1): 283-292.
*/

例子:

//这里以小木虫的{三维曲面绘制求助<http://emuch.net/bbs/viewthread.php?tid=5999724&fpage=1>}例子为例,下面的data即来自帖子,依次回车执行各条命令可看到结果


x=data(,1);//获得data第1列数据
y=data(,2);
z=data(,3);
a=max(x);
b=min(x);
x0=MeshGrid2(b,a,50);//生成x坐标网格
a=max(y);
b=min(y);
y0=MeshGrid2(b,a,50);//生成y坐标网格
y0=y0';//注意要转置一次
z0=InterpolationRBF2D(x,y,z,x0,y0);//调用本函数求解
cloud(x0,y0,z0)//可以绘制云图查看效果
{<http://img2.ph.126.net/BszLgWAey2XLqMujFX_cNg==/1059753287333311586.png>}
