多项式插值(GNUOctave(版本9.4.0))

28.5多项式插值¶

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

: 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可以直接与一起使用polyval使用拟合的多项式来估计值。

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

‘yf’: 的每个值的多项式值x.

第二输出可从polyval以计算预测值的统计误差极限。特别是的标准偏差p系数从

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

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

f = @(x) x.^2 + 5; # data-generating function x = 0:5; y = f (x); ## Fit data to polynomial A*x^3 + B*x^1 p = polyfit (x, y, logical ([1, 0, 1, 0])) ⇒ p = [ 0.0680, 0, 4.2444, 0 ] ## Fit data to polynomial using all terms up to x^3 p = polyfit (x, y, 3) ⇒ p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ]

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

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

x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
## Uniform breaks
breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
pp1 = splinefit (x, y, breaks);
## Breaks interpolated from data
pp2 = splinefit (x, y, 10);  % 11 breaks, 10 pieces
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})

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

## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
## Piecewise constant
pp1 = splinefit (x, y, 8, "order", 0);
## Piecewise linear
pp2 = splinefit (x, y, 8, "order", 1);
## Piecewise quadratic
pp3 = splinefit (x, y, 8, "order", 2);
## Piecewise cubic
pp4 = splinefit (x, y, 8, "order", 3);
## Piecewise quartic
pp5 = splinefit (x, y, 8, "order", 4);
## Plot
xx = linspace (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])
axis tight
ylim auto
legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})

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

## Data (100 points)
x = 2 * pi * [0, (rand (1, 98)), 1];
y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
## No constraints
pp1 = splinefit (x, y, 10, "order", 5);
## Periodic boundaries
pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "no constraints", "periodic"})

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

## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (2 * x) + 0.1 * randn (size (x));
## Breaks
breaks = linspace (0, 2 * pi, 10);
## Clamped endpoints, y = y' = 0
xc = [0, 0, 2*pi, 2*pi];
cc = [(eye (2)), (eye (2))];
con = struct ("xc", xc, "cc", cc);
pp1 = splinefit (x, y, breaks, "constraints", con);
## Hinged periodic endpoints, y = 0
con = struct ("xc", 0);
pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "clamped", "hinged periodic"})

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

## Data
x = linspace (0, 2*pi, 200);
y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
## Add outliers
x = [x, linspace(0,2*pi,60)];
y = [y, -ones(1,60)];
## Fit splines with hinged conditions
con = struct ("xc", [0, 2*pi]);
## Robust fitting, beta = 0.25
pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
## Robust fitting, beta = 0.75
pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
## No robust fitting
pp3 = splinefit (x, y, 8, "constraints", con);
## Plot
xx = linspace (0, 2*pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
plot (x, y, ".", xx, [y1; y2; y3])
legend ({"data with outliers","robust, beta = 0.25", ...
         "robust, beta = 0.75", "no robust fitting"})
axis tight
ylim auto

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

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合并为一个函数。这些函数中的每一个都用邻接区间表示。