基于MATLAB的B样条插值拟合算法与分段多项式（附完整代码） |
您所在的位置：网站首页 › 拟合图像怎么做 › 基于MATLAB的B样条插值拟合算法与分段多项式（附完整代码） |
基于MATLAB的B样条插值拟合算法与分段多项式（附完整代码）
|
一. B样条函数
B样条函数的MATLAB代码如下： S=spapi(k,x,y) %k为用户选定的B样条阶次，一般以4和5居多 例题1分别用B样条函数对y和f(x)中的自选数据进行5次B样条函数拟合，并与三次分段多项式样条函数拟合的结果相比较。 解： MATLAB代码如下： clc;clear; %%y函数部分 x0=[0,0.4,1,2,pi]; y0=sin(x0); ezplot('sin(t)',[0,pi]); hold on %三次分段多项式样条插值 sp1=csapi(x0,y0); fnplt(sp1,'--'); %5次B样条插值 sp2=spapi(5,x0,y0); fnplt(sp2,':') %%f(x)函数部分 x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); figure; ezplot('(x.^2-3*x+5).*exp(-5*x).*sin(x)',[0,1]), hold on sp1=csapi(x,y); fnplt(sp1,'--'); sp2=spapi(5,x,y); fnplt(sp2,':')运行结果：
在MATLAB中，可以利用以下函数来求S的k阶导数： Sd=fnder(S,k)如果是求多变量函数的偏导数，则格式如下： Sd=fnder(S,[k1,...,kn]) 例题2任取函数f(x)的一些数据点，分别用三次分段多项式样条函数与B样条插值函数进行拟合，求出该函数的导数，并与理论推导结果进行比较。 解： MATLAB代码如下： clc;clear; syms x; f=(x^2-3*x+5)*exp(-5*x)*sin(x); ezplot(diff(f),[0,1]) %理论结果 hold on, x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); sp1=csapi(x,y); %建立三次样条函数 dsp1=fnder(sp1,1); fnplt(dsp1,'--'); %绘制样条图 hold on, sp2=spapi(5,x,y); %5阶次B样条插值 dsp2=fnder(sp2,1); fnplt(dsp2,':'); axis([0,1,-0.8,5])运行结果：
自行选择z=f(x,y)中的数据，来拟合 解： MATLAB代码如下： clc;clear; %拟合曲面 x0=-3:.3:3; y0=-2:.2:2; [x,y]=ndgrid(x0,y0); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); sp=spapi({5,5},{x0,y0},z); %B样条 dspxy=fnder(sp,[1,1]); fnplt(dspxy) %生成样条图 %理论方法 syms X Y; Z=(X^2-2*X)*exp(-X^2-Y^2-X*Y); figure; ezsurf(diff(diff(Z,X),Y),[-3 3],[-2 2]) %对符号变量表达式做三维表面图运行结果：
S作为样条函数，对其进行数值积分格式如下： f=fnint(S) 例题4考虑以下积分中较稀疏的样本点，用样条积分的方式求出定积分及积分函数。 解： MATLAB函数如下： clc;clear; x=[0,0.4,1,2,pi]; y=sin(x); %建立三次样条函数并积分 sp1=csapi(x,y); a=fnint(sp1,1); xx1=fnval(a,[0,pi]); integral1=xx1(2)-xx1(1) %建立B样条函数并积分 sp2=spapi(5,x,y); b=fnint(sp2,1); xx2=fnval(b,[0,pi]); integral2=xx2(2)-xx2(1) %绘制曲线 ezplot('-cos(t)+2',[0,pi]); %不定积分可以上下平移 hold on, fnplt(a,'--'); fnplt(b,':');运行结果：
根据间断数和系数生成分段多项式pp，每个系数的值都是长度为d的向量。MATLAB格式如下： pp=mkpp(breaks,coefs,d)可以使用ppval来计算特定点处的分段多项式，也可以使用unmkpp来提取有关分段多项式的详细信息。 例题5创建任意一个分段多项式，使得它在区间[0,4]内具有三次多项式，在区间[4,10]内具有二次多项式，在区间[10,15]内具有四次多项式。 解： MATLAB代码如下： clc;clear; breaks=[0 4 10 15]; coefs=[0 1 -1 1 1;0 0 1 -2 53;-1 6 1 4 77]; %一共有三段 pp=mkpp(breaks,coefs) %画图 xq=0:0.01:15; plot(xq,ppval(pp,xq)) line([4 4],ylim,'LineStyle','--','Color','k') line([10 10],ylim,'LineStyle','--','Color','k')运行结果：
创建一个具有四个区间的单个分段多项式，这些区间在两个二次多项式之间交替。 解： MATLAB代码如下： clc;clear; %显示一个二次多项式在[-8,-4]区间上的结果 subplot(2,2,1) cc=[-1/4 1 0]; pp1=mkpp([-8 -4],cc); xx1=-8:0.1:-4; plot(xx1,ppval(pp1,xx1),'k-') %在[-4,0]区间上的求反 subplot(2,2,2) pp2=mkpp([-4 0],-cc); xx2=-4:0.1:0; plot(xx2,ppval(pp2,xx2),'k-') %将二次多项式扩展到四个区间形成的分段多项式 %显示一阶导数，该导数利用unmkpp分解分段多项式构造而成 subplot(2,1,2) pp=mkpp([-8 -4 0 4 8],[cc;-cc;cc;-cc]); xx=-8:0.1:8; plot(xx,ppval(pp,xx),'k-') [breaks,coefs,k,d]=unmkpp(pp); dpp=mkpp(breaks,repmat(k-1:-1:1,d*1,1).*coefs(:,1:k-1),d);运行结果：
|
【本文地址】