matlab样条插值拟合,科学网—样条函数插值拟合 - 李继存的博文
样条函数插值拟合
2014–02–11 09:26:49
在拟合势能函数的时候, 除解析式外, 也可以利用样条函数进行拟合. 样条拟合与其插值正好相反: 已知函数在节点上的值求任意位置的值, 做插值; 已知函数的某些组合值求函数节点上的值, 属于拟合. 由于样条函数可以化为节点函数值的线性表达式, 这样就可以将待求参数线性化, 得到最优情况下函数的形状, 为寻找合适的解析式提供依据, 当然也可以直接利用得到的离散数据拟合解析式.
样条函数可以是零阶, 一阶, 二阶, 三阶或更高阶. 实际使用中, 三阶使用最为普遍. 由于三次样条的构造需要求解一个三对角线性方程组, 其显式解很难得到, 所以线性化结果很繁琐.零阶
在每一区间上样条函数为常量, 函数整体呈台阶状. 对等距情况, 计算时最好使用就近原则, 取最近点的值作为拟合点, 可用i=nint(x/dx)实现.
一阶
在每一区间上样条函数为线性函数, 函数整体呈折线状$.$
fi(x)=yi+yi+1−yiΔx(x−xi)
此式自动满足函数值连续条件, 即零阶连续.
二阶
在每一区间上, 样条函数为二次函数, 整体一阶连续, 即有连续的导数. 但仅有节点的函数值不能唯一确定整个函数, 还须提供某一节点上的导数值, 一般可令端点的导数值为零. 二次样条函数在偶数点的曲率不连续. 由二阶开始, 插值函数不再具有局域性, 改变某一节点, 函数整体都会改变. 使得线性系数分离很困难.
fi(x)ki+1=yi+ki(x−xi)+ki+1−ki2Δx(x−xi)2=−ki+2yi+1−yiΔx,ki=2yi+1−yiΔx−ki+1
可化简得
fi(x)α=yi+ki(x−xi)+(yi+1−yi−kiΔx)(x−xiΔx)2=α2yi+1+(1−α2)yi+(1−α)ωki=ω/Δx,ω=x−xi
对势能函数, 一般满足远距离处导数为零, 故可使用自然条件 kn=0 , 由此, 可推知所有系数 ki .
令 ki=2ΔxTi,Δi=yi+1−yi, 则 Ti 满足递推式
Ti=Δi−Ti+1
可求得
Ti=∑j=in−1(−1)j−iΔj
样条函数可写为
fi(x)=α2yi+1+(1−α2)yi+2α(1−α)Ti,α=(x−xi)/Δx
对不等距划分, 令 Δi=yi+1−yixi+1−xi, ki 满足如下递推式
ki=2Δi−ki+1
求得
ki=2∑j=in−1(−1)j−iΔj
三阶
对等距划分的均匀样条, 设节点为 1,2,....n, 若 x∈[xi,xi+1],a=x−xi,b=xi+1−x,a+b=h=Δx, 则
6hfi(x)=6(ayi+1+byi)+a(a2−h2)Mi+1+b(b2−h2)Mi
Mi 为节点的二阶导数, 对应于力学上的弯矩, 满足下面的方程
Mi+4Mi+1+Mi+2=di+1=6h2(yi+2−2yi+1+yi),i=1,2...n−2
要求的 Mi 个数为 n, 而对应的方程数目为 n−2, 故还需两个边界条件才能唯一确定, 边界条件可取为两端点的导数值或是二阶导数值. 常用的自然边界条件指 M1=Mn=0. 加上边界条件后便可求得 M2,M3,...Mn−1
Mi 满足的方程为
0M2M3⋮Mn−3Mn−2+++⋮++4M24M34M4⋮4Mn−24Mn−1+++⋮++M3M4M5⋮Mn−10===⋮==d2d3d4dn−2dn−1
写为矩阵形式
⎡⎣⎢⎢⎢⎢⎢41⋮……14⋮1001⋮41……⋮14⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢M2M3⋮Mn−2Mn−1⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢d2d3⋮dn−2dn−1⎤⎦⎥⎥⎥⎥⎥⎥
可见, 此方程为三对角方程组, 对角占优, 存在唯一解, 可利用所谓的追赶法求解. 中文常称的追赶法, 是Thomas方法的形象翻译. 大致求解过程分为两步:追: 利用消元法将原方程化为二对角方程, 向前递推, 使系数矩阵主对角线变为1. 由第一个方程, 得到 M2 和 M3 的关系, 将其带入第二个方程, 消去主对角线下方系数. 以此进行, 最终追到 Mn−1, 得到其解. 变换后, 其方程为
⎡⎣⎢⎢⎢⎢⎢⎢10⋮……A21⋮000A3⋮10……⋮An−21⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢M2M3⋮Mn−2Mn−1⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢D2D3⋮Dn−2Dn−1⎤⎦⎥⎥⎥⎥⎥⎥
赶: Mn−1 已经追得, 然后由此倒推, 得到其他 Mi 值.
对上面的方程, 由于系数是固定的, Thomas方法的递推式为
A2D2Mn−1=14,=A2d2,=Dn−1,AiDiMi=14−Ai−1=Ai(di−Di−1)=Di−AiMi+1
对 Ai, 可求得其通式
Ai=2−3√ti+1ti−1,t=(2+3√)2
对 Di, 向后递推至 D2, 一般项可写为
Dj=∑k=2j(−1)j−kPjkdk
对 Mi, 向前递推至 Mn−1, 一般项可写为
Mi=∑j=in−1(−1)j−iPj−1iDj
综合上面两个结果, 得到
MiPnm=∑j=in−1∑k=2j(−1)2j−i−kPj−1iPjkdk=∑j=in−1∑k=2j(−1)i+kPj−1iPjkdk,i=2,3,...n−1=∏l=mnAl=AmAm+1…An−1An,n>m
根据插值公式
dk6hfi(x)=6h2(yk+1−2yk+yk−1)=6(αyi+1+βyi)+α(α2−h2)Mi+1+β(β2−h2)Mi
以划分间距 h 为单位, 约化上述公式
fi(x)α=αyi+1+βyi+α(α2−1)μi+1+β(β2−1)μi=ah,β=bh=1−α,μi=Mi6/h2
由此可见, 虽然三次样条函数仍可写为节点值的线性形式, 但其系数十分复杂.
上面公式看起来清楚, 但是实际计算时需要计算 Mi 和 Mi+1, 两个三重循环, 整体计算量为 O(N3). 利用 Ai 的近似关系和 Mi 的递推关系可以将公式的计算量减少一些.
fi(x)=αyi+1+βyi+α(α2−1)μi+1+β(β2−1)μi=αyi+1+βyi+β(β2−1)Di+[α(α2−1)−β(β2−1)Ai]μi+1
对不等距划分, 上面的递推公式太过复杂, 很难写出一般项了. 样条函数可写为
fi(x)hi=ayi+1+byihi+a(a2−h2)Mi+1+b(b2−h2)Mi6hi=xi+1−xi,a=x−xi,b=hi−a
Mi 满足方程
hiMi+2(hi+hi+1)Mi+1+hi+1Mi+2=6(yi+2−yi+1hi+1−yi+1−yihi),i=1,2,…,n−2
代码及测试测试结果
awk实现的代码如下awk ' BEGIN{ Ndat=0 }
NF==3 { Ndat++; X[Ndat]=$2; Y[Ndat]=$3 }
END {
for(i=1; i
for(i=2; i
A[1]=0; A[Ndat]=0
D[1]=0; D[Ndat]=0
for(i=2; i<=Ndat-1; i++) {
ai=dX[i-1]
bi=2*(dX[i-1]+dX[i])
ci=dX[i]
A[i]=ci/(bi-ai*A[i-1])
D[i]=(d[i]-ai*D[i-1])*A[i]/ci
}
M[Ndat]=0
for(i=Ndat-1; i>1; i--) M[i]=D[i]-A[i]*M[i+1]
# for(i=1; i<=Ndat; i++) print i, dX[i], d[i], A[i], D[i], M[i]
h=X[2]-X[1]
for(x=X[1]; x
}
function SP2(x, i, Ndat, X, Y, dX, j,a,ki) {
if(i==0) {
for(i=1; ix) break
}
ki=0
for(j=i; j<=Ndat-1; j++) ki += 2*(-1)^(j-i)*(Y[j+1]-Y[j])/dX[j]
a=(x-X[i])/dX[i]
return a^2*Y[i+1]+(1-a^2)*Y[i]+(1-a)*ki*(x-X[i])
}
#
function SP3(x, i, Ndat, X, Y, dX, M, a,b,h) {
if(i==0) {
for(i=1; ix) break
}
h=dX[i]
a=x-X[i]; b=h-a
return (a*Y[i+1]+b*Y[i])/h + ( a*(a^2-h^2)*M[i+1]+b*(b^2-h^2)*M[i] )/(6*h)
}
#
function Psp3(Ndat, i,j,k,a) {
a=(2+sqrt(3))^2
for(i=2; i
for(j=i; j
P[i,j]=1
for(k=i; k<=j; k++) P[i,j] *= 2-sqrt(3)*(a^k+1)/(a^k-1)
}
}
}
#
function uniSP3(x, i, Ndat, X, Y, j,k,a,b,h) {
if(i==0) {
for(i=1; ix) break
}
h=X[2]-X[1]
a=(x-X[i])/h; b=1-a
for(j=1; j<=Ndat; j++) Coef[j]=0
Coef[i] += b
Coef[i+1] += a
Rsec=b*(b^2-1)
for(k=2; k<=i; k++) {
Rtmp=Rsec * (-1)^(i-k) * P[k,i]
Coef[k+1] += Rtmp
Coef[k] += -2*Rtmp
Coef[k-1] += Rtmp
}
Rsec=a*(a^2-1)-b*(b^2-1)*P[i,i]
for(j=i+1; j<=Ndat-1; j++) {
Rij=Rsec * (-1)^(i+1)
if(j!=i+1) Rij *= P[i+1,j-1]
for(k=2; k<=j; k++) {
Rtmp=Rij * (-1)^k * P[k,j]
Coef[k+1] += Rtmp
Coef[k] += -2*Rtmp
Coef[k-1] += Rtmp
}
}
F=0
for(j=1; j<=Ndat; j++) F += Coef[j]*Y[j]
return F #( a*Y[i+1]+b*Y[i] + a*(a^2-1)*Mi(i+1, Ndat) + b*(b^2-1)*Mi(i, Ndat) )
}
#
function Mi(i, Ndat, j,k,Rtmp,ret) {
ret=0
for(j=i; j<=Ndat-1; j++) {
Rtmp=(-1)^i
if(j!=i) Rtmp *= P[i,j-1]
for(k=2; k<=j; k++)
ret += Rtmp*(-1)^k*P[k,j] *d[k]
}
return ret
}
' ABS >ABS.sp
# ' CUB > CUB.sp
#' LJ >LJ.sp
利用Matlab的csape函数可以进行三次样条函数的插值, 示例代码如下format long
x=1:0.02:15;
y=-1E3*(-12./x.^13+6./x.^7);
pp=csape(x,y,'second',[0,0]);
pp=csape(x,y);
xsp=1:0.01:1.5;
ysp=ppval(pp,xsp);
yst=-1E3*(-12./xsp.^13+6./xsp.^7);
plot(x,y,'o',xsp,ysp,'-',xsp,yst,'-')
axis([1 1.5 -600 600])
FID=fopen('LJ.mat', 'w');
Ndat=length(xsp);
for i=1:Ndat
fprintf(FID, '%f %fn', xsp(i), ysp(i));
end
几个测试函数的结果
样条函数的积分
对已经拟合的样条函数进行数值积分时可利用可利用Simpson方法. Simpson方法具有三阶精度, 对不超过三次的多项式精确成立, 很适合于二次和三次样条函数.
转载本文请联系原作者获取授权,同时请注明本文来自李继存科学网博客。
链接地址:http://blog.sciencenet.cn/blog-548663-775279.html
上一篇:GAMESS编译使用简记
下一篇:Bash脚本中使用颜色
matlab样条插值拟合,科学网—样条函数插值拟合 - 李继存的博文相关推荐
- java 处理pdb文件格式_科学网—PDB文件格式说明 - 李继存的博文
2015-06-05 20:31:19 2017-01-22 20:09:21 据参考资料增补 PDB(Protein Data Bank)是一种标准文件格式, 其中包含原子的坐标等信息, 提交给 P ...
- markdown 流程图js_科学网—让Markdown支持ASCII流程图和JavaScript流程图 - 李继存的博文...
2014-12-25 12:08:34 计算机领域中一直存在两种不同的理念并彼此竞争, 可视化与可控化, 或称为所见即所得与所愿即所得. 前者是Windows的典型做法, 而后者是Linux的典型理念 ...
- linux子系统安装gromacs,科学网—Windows下GROMACS程序的编译 - 李继存的博文
2015-12-07 22:12:05 总的来说, Windows下的GROMACS程序用于模拟意义不大, 对于长时间的模拟, 我都是放在Linux服务器上进行的. 但将Windows下的GROMAC ...
- 基于MATLAB的二维与三维插值拟合运算(附完整代码)
· 一. 一维插值 interp1函数在上个博客中(如下链接)已经更新了,此处再补充两个相关例题. 基于MATLAB的数据插值运算:Lagrange与Hermite算法(附完整代码)_唠嗑!的博客-C ...
- matlab基本矩阵运算,科学网—matlab中矩阵基本运算 - 成爱芳的博文
以矩阵A为例: (1)转置矩阵求取 AT transpose(A) >> A=[1 0 3; 2:4; 3 1 0] A = 1 0 3 2 3 4 3 ...
- matlab画复变函数,科学网—复数复变函数的Matlab计算与绘图 - 周铁戈的博文
复数复变函数的Matlab计算与绘图 周铁戈 复数的表示 存在两种表示方法,一种是代数式,一种是指数式,在Matlab中的方式如下: >> z=1+2i #代数式,1 ...
- matlab聚类算法,科学网—matlab-聚类算法笔记 - 孙月芳的博文
MATLAB提供了两种方法进行聚类分析: 1.利用clusterdata 函数对数据样本进行一次聚类,这个方法简洁方便,其特点是使用范围较窄,不能由用户根据自身需要来设定参数,更改距离计算方法: 2. ...
- matlab 水平投影,科学网—Matlab中如何将投影信息写入到shape文件中 - 朱永超的博文...
在Matlab中保存shape格式数据时,没有具体的函数可以将投影信息直接写入到shape文件中,不过可以通过另外一种方式实现.看下shape格式的文件不难发现,shape文件的投影信息是一个单独的文 ...
- matlab流量结构分析,科学网-分享求解“结构分解分析(SDA)”各项均值的MATLAB程序-计军平的博文...
点此下载(MATLAB File Exchange) [2015.02.18补充]其他研究人员的MATLAB代码 Sayago-Gomez, Juan Tomas, (2014), Matlab Co ...
最新文章
- 目录页码错误未定义书签怎么解决_目录页码对不齐应该怎么办?这2种方法,工作效率大增...
- 揭秘一线互联网企业 前端JavaScript高级面试
- PyTorch框架学习七——自定义transforms方法
- 笨方法“学习python笔记之列表
- 阅读《大型网站技术架构》前两章心得体会及总结
- nginx: [emerg] bind() to 0.0.0.0:66 failed (98: Address already in use)
- matlab 神经网络工具箱的实用
- AD09之与AD6版本使用不同对比
- Android Gradle构建脚本
- 地球的3D模型制作教程【3DsMax】
- ThinkBook 14 G2 ITL 重装系统 笔记
- Linux之RPM包的命名规则和包的依赖性
- vue引入 wps在线编辑版,可进行 预览,编辑, 打印等功能。
- 最后一批90后开始养生了,中医科普短视频会火吗?
- BoardCast广播组件
- ffmpeg 图片序列转视频
- 计算机操作系统-文件管理 知识点归纳
- VS2017环境下GMap的学习及开发(一)
- 基于Opencv的颜色识别
- 考考眼力,找下代码的bug