28.5多项式插值¶

Octave为各种插值提供了良好的支持,其中大部分在中进行了描述插值上述章节中描述的函数的一个简单替代方案是将单个多项式或分段多项式(样条曲线)拟合到一些给定的数据点。为了避免高度波动的多项式,人们通常希望将低阶多项式拟合到数据中。这通常意味着在最小二乘意义上拟合多项式是必要的,这就是拟合函数确实如此。

: p = polyfit (x, y, n) ¶

: [p, s] = polyfit (x, y, n) ¶

: [p, s, mu] = polyfit (x, y, n) ¶

返回多项式的系数</p>(x)的n使点拟合的最小二乘误差最小化[x(:), y(:)].

n通常是一个≥0的整数,用于指定近似多项式的阶数。如果n是一个逻辑向量,它被用作amask,以选择性地强制使用或忽略相应的多项式系数。

多项式系数在行向量中返回</p>.输出</p>可以直接与一起使用多项式求值使用拟合的多项式来估计值。

可选输出s是一个包含以下字段的结构体:

‘yf’

的每个值的多项式值x.

‘十、’

用于计算多项式效率的范德蒙德矩阵。

‘R’

来自QR分解的三角因子R。

‘C’

未缩放的协方差矩阵,形式上等于的倒数x*x,但是以最小化舍入误差传播的方式进行计算。

‘df’

自从度。

‘normr’

残差的范数。

第二输出可从多项式求值以计算预测值的统计误差极限。特别是的标准偏差</p>系数从

sqrt(diag(s.C)/s.df) * s.normr.

当第三输出,μ,对原始数据进行了居中和缩放,这可以提高拟合的数值稳定性。系数</p>与中的多项式相关联

xhat= (x- μ1.μ2.
这里的μ(1) =平均值(x和μ(2) =标准(x).

示例1:逻辑n和整数n

f=@(x)x.^2+5;#数据生成函数x=0:5;y=f(x);##将数据拟合到多项式A*x^3+B*x^1p=多元拟合(x,y,逻辑([1,0,1,0]))⇒ p=[0.0680,0,4.2444,0]##使用所有项将数据拟合到多项式,直到x^3p=polyfit(x,y,3)⇒ p=[4.9608e-17,1.0000e+00,-3.3813e-15,5.0000e+00]

详见: 多项式求值, 聚烷烃, 根, 范德, zscore.

在单个多项式不够好的情况下,解决方案是使用几个多项式拼凑在一起。函数花键配合将分段多项式(样条曲线)拟合到一组数据。

: pp = splinefit (x, y, breaks) ¶

: pp = splinefit (x, y, p) ¶

: pp = splinefit (…, "periodic", periodic) ¶

: pp = splinefit (…, "robust", robust) ¶

: pp = splinefit (…, "beta", beta) ¶

: pp = splinefit (…, "order", order) ¶

: pp = splinefit (…, "constraints", constraints) ¶

拟合带断点(节)的分段三次样条曲线中断到那时的noisy数据,x和y.

x是一个向量,并且y是向量或N-D数组。如果y是一个N-D数组,那么x(j) 与匹配yj

</p>是一个正整数,定义沿x和</p>+1是中断次数。每个区间的点数相差不超过1。

可选属性周期性是一个逻辑值,用于指定周期边界条件是否应用于样条曲线。周期的长度为最大(中断)-分钟(中断)。默认值为错误的.

可选属性强健的是一个逻辑值,用于指定是否应用稳健拟合来减少输出数据点的影响。执行三次加权最小二乘迭代。权重是根据前一个残差计算得出的。外部识别的灵敏度从属性控制贝塔。的值贝塔被限制在范围内,0<贝塔1.默认值为贝塔= 1/2. 接近0的值赋予所有数据相等的权重。的价值不断增加贝塔减少数据泄露的影响。接近1的值可能会导致不稳定或秩效率。

拟合的样条作为分段多项式返回,pp,可以使用进行评估ppval.

样条曲线是从次数为的多项式构成的顺序默认值为立方体,顺序=3.具有P个片的样条具有P+顺序自从度。在周期性边界条件下,自从度降为P。

可选属性,约束,是一个在拟合上指定线性约束的结构体。该结构体具有三个字段,xc,yc和抄送.

xc

约束的x位置的向量。

yc

约束位置处的值xc。默认值是一个零数组。

抄送

系数(矩阵)。默认值是一个一的数组。行的数量被限制为分段多项式的阶数,顺序.

约束是0到的阶导数的线性组合顺序-1根据

cc(1,j)*y(xc(j))+cc(2,j)*y'(xc)+…=yc(:,…,:,j)。

详见: interp1, unmkpp, ppval, 样条曲线, pchip, ppder, ppint, ppjumps.

的数量中断(或节点)是抑制输入数据中存在的噪声的重要因素,x和y。下面的例子说明了这一点。

x=2*pi*rand(1200);y=sin(x)+sin(2*x)+0.2*randn(大小(x));##均匀破裂=林空间(0,2*pi,41);%41个断裂,40个spp1=花键配合(x,y,断裂);##根据数据插值的断点pp2=样条曲线拟合(x,y,10);%11个断点,10个##Plotxx=林空间(0,2*pi,400);y1=ppval(pp1,xx);y2=ppval(pp2,xx);plot(x,y,“.”,xx,[y1;y2])axis tightylim autolegend({“data”,“41个断裂,40个”,“11个断裂,10个”})

其结果可以在中看到图28.1.

图28.1:具有41个断点的分段多项式与具有11个断点的多项式的拟合的比较。具有大量断裂的拟合表现出快速的波纹,这在基础函数中是不存在的。

分段多项式拟合,从提供花键配合,具有连续的导数,最高可达顺序-例如,一个三次fit具有连续的一阶导数和二阶导数。代码对此进行了演示

##数据(200点)x=2*pi*rand(1200);y=sin(x)+sin(2*x)+0.1*randn(大小(x));##分段常数pp1=花键配合(x,y,8,“顺序”,0);##分段线性pp2=花键配合(x,y,8,“顺序”,1);##分段象限pp3=花键配合(x,y,8,“顺序”,2);##分段cubiccp4=花键配合(x,y,8,“order”,3);##分段四次pp5=花键配合(x,y,8,“顺序”,4);##Plotxx=林空间(0,2*pi,400);y1=ppval(pp1,xx);y2=ppval(pp2,xx);y3=ppval(pp3,xx);y4=ppval(pp4,xx);y5=ppval(pp5,xx);plot(x,y,“.”,xx,[y1;y2;y3;y4;y5])轴紧密型自动egend({“data”,“order 0”,“order 1”,“order2”,“Order3”,“order 4”})

其结果可以在中看到图28.2.

图28.2:具有8个中断的分段常数、线性、二次、三次和四次多项式与噪声数据的比较。高阶解更准确地表示底层函数,但也会带来计算复杂性的代价。

当提供拟合的基础函数是周期性的时,花键配合能够应用表现周期性拟合所需的边界条件。下面的代码演示了这一点。

##数据(100点)x=2*pi*[0,(rand(1,98)),1];y=sin(x)-cos(2*x)+0.2*randn(大小(x));##无约束spp1=花键配合(x,y,10,“顺序”,5);##周期性边界esp2=花键配合(x,y,10,“顺序”,5,“周期性”,true);##Plotxx=林空间(0,2*pi,400);y1=ppval(pp1,xx);y2=ppval(pp2,xx);plot(x,y,“.”,xx,[y1;y2])轴紧密型自动生成({“数据”,“无约束”,“周期性”})

其结果可以在中看到图28.3.

图28.3:分段多项式在有和无周期边界条件的情况下适用于有噪声的周期函数的比较。

也可以添加更复杂的约束。例如,下面的代码将使用在端点处夹紧的值来计算一个周期性拟合,以及在端点处铰接的第二个周期性配合。

##数据(200点)x=2*pi*rand(1200);y=sin(2*x)+0.1*randn(大小(x));##Breaksbreaks=林空间(0,2*pi,10);##箝位端点,y=y'=0xc=[0,0,2*pi,2*pi];cc=[(eye(2)),(eye(1))];con=结构体(“xc”,xc,“cc”,cc);pp1=花键配合(x,y,断裂,“约束”,con);##铰接的周期性端点,y=0con=结构体(“xc”,0);pp2=花键配合(x,y,breaks,“constraints”,con,“periodical”,true);##Plotxx=林空间(0,2*pi,400);y1=ppval(pp1,xx);y2=ppval(pp2,xx);plot(x,y,“.”,xx,[y1;y2])轴紧密型自动校准({“数据”,“夹紧”,“铰接周期”})

其结果可以在中看到图28.4.

图28.4:两个周期分段三次拟合与噪声周期信号的比较。一个配合的端点被夹住,第二个配合的终点被磨光。

这个花键配合函数还提供了强健的拟合,其中减少了外围数据的影响。在下面的示例中,提供了三种不同的配合。两个具有不同的外部抑制水平,第三个说明了非鲁棒解决方案。

##Datax=林空间(0,2*pi,200);y=sin(x)+sin(2*x)+0.05*randn(大小(x));##添加outliersx=[x,林空间(0,2*pi,60)];y=[y,-ones(1.60)];##使用铰接条件拟合样条曲线scon=struct(“xc”,[0,2*pi]);##稳健拟合,β=0.25pp1=花键拟合(x,y,8,“约束”,con,“β”,0.25);##稳健拟合,β=0.75pp2=花键拟合(x,y,8,“约束”,con,“β”,0.75);##无稳健配合pp3=花键配合(x,y,8,“约束”,con);##Plotxx=林空间(0,2*pi,400);y1=ppval(pp1,xx);y2=ppval(pp2,xx);y3=ppval(pp3,xx);plot(x,y,“.”,xx,[y1;y2;y3])图例({“带异常值的数据”,“稳健,β=0.25”,…“稳健,贝塔=0.75”,“无稳健拟合”})轴紧线

其结果可以在中看到图28.5.

图28.5:两种不同稳健拟合水平的比较(贝塔=0.25和0.75)到与外围数据组合的噪声数据。传统合身,不合身(贝塔=0)。

一种非常具体的多项式解释形式是帕德近似。对于控制系统,连续时间延迟可以用近似值非常简单地建模。

: [num, den] = padecoef (T) ¶

: [num, den] = padecoef (T, N) ¶

计算N连续时滞的一阶Padé逼近T以传递函数的形式。

的Padé近似exp(-sT)从以下方程定义

Pn(s)exp(-sT)~--------Qn(s)

其中Pn(s)和Qn(s)都是N从以下表达式定义的次有理函数

N(2N-k)！NkPn(s)=总和---------------(-sT)k=0(2N)！k(N-k)！Qn(s)=Pn(-s)

输入T和N必须是非负数字标量。如果N未指定,则默认为1。

输出行向量号码和兽穴包含以s的降序排列的分子和分母系数。两者都是N阶多项式。

例如

t=0.1;n=4;[num,den]=padecoef(t,n)⇒ num=1.0000e-04-2.0000e-02 1.8000e+00-8.4000e+01 1.6800e+03⇒ den=1.0000e-04 2.0000e-02 1.8000e+00 8.4000e+01 1.6800e+03

函数,ppval,计算分段多项式,从创建mkpp或其他方式,以及unmkpp返回有关分段多项式的详细信息。

以下示例显示如何将两个线性函数和aquadratic合并为一个函数。这些函数中的每一个都用邻接区间表示。

x=[-2,-1,1,2];p=[0,1,0;1,-2,1;0,-1,1];pp=mkpp(x,p);xi=林空间(-2,2,50);yi=ppval(pp,xi);情节(xi、易);

: pp = mkpp (breaks, coefs) ¶

: pp = mkpp (breaks, coefs, d) ¶

从采样点构造分段多项式(pp)结构体中断和系数系数.

中断必须是值严格递增的向量。intervals的数量从镍长中断1..

什么时候m是多项式阶系数大小必须为:镍通过m1..

的第i行系数, 系数(我,:),包含多项式在我-第th个间隔,按最高顺序排列(m)至最低(0度

系数也可以是多维数组,指定avector值或数组值多项式。在这种情况下,多项式阶m从的最后一个维度的长度定义系数。第一个维度的大小从标量或向量给出d如果d未给定它设置为1.在这种情况下</p>(r, 我, :)包含的系数r-区间上定义的第th个多项式我.无论如何系数被重塑为大小为2D的矩阵[镍戳d) m].

编程说明:ppval在处计算多项式xi- 中断(i),即从中减去当前间隔的下端xi。在创建具有mkpp.

详见: unmkpp, ppval, 样条曲线, pchip, ppder, ppint, ppjumps.

: [x, p, n, k, d] = unmkpp (pp) ¶

提取分段多项式结构体的分量pp.

此函数与的相反mkpp:它提取输入到mkpp需要创建分段多项式结构体pp。下面的代码明确了这种关系:

[breaks,coeffes,numinter,order,dim]=unmkpp(pp);pp2=mkpp(断裂、系数、dim);

分段多项式结构体第2页以这种方式获得,与原件不一样pp.同样可以通过直接访问结构体的字段来获得pp.

组件包括:

x

采样点或断点。

</p>

采样区间中点的多项式系数。</p>(我, :)包含多项式上区间的系数我从最高到最低排序。如果d1.然后</p>是一个大小矩阵[n戳d) m],其中我1.d)行是所有的系数d区间中的多项式我.

n

多项式片段或区间的数量,n长x1..

k

多项式的阶数加1。

d

为每个区间定义的多项式数。

详见: mkpp, ppval, 样条曲线, pchip.

: yi = ppval (pp, xi) ¶

评估分段多项式结构体pp在点上xi.

如果pp描述了一个标量多项式函数,结果是一个与相同形状的数组xi。否则,结果的大小为[pp.dim,长度(xi)]如果xi是向量,或者[pp.dim,尺寸(xi)]如果它是多维数组。

详见: mkpp, unmkpp, 样条曲线, pchip.

: ppd = ppder (pp) ¶

: ppd = ppder (pp, m) ¶

计算分段m-分段多项式结构体的th导数pp.

如果m则计算一阶导数。

详见: mkpp, ppval, ppint.

: ppi = ppint (pp) ¶

: ppi = ppint (pp, c) ¶

计算分段多项式结构体的积分pp.

c,如果给定,则为积分常数。

详见: mkpp, ppval, ppder.

: jumps = ppjumps (pp) ¶

评估分段多项式的边界跳跃。

如果有n区间,以及的维度pp是d,得到的数组具有维度[d,n-1].

详见: mkpp.