基于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])

运行结果：

例题3

自行选择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')

运行结果：

例题6

创建一个具有四个区间的单个分段多项式，这些区间在两个二次多项式之间交替。

解：

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);

 运行结果：