octave: Polynomial Interpolation
28.5 Polynomial Interpolation
=============================
Octave comes with good support for various kinds of interpolation, most
of which are described in 
Interpolation. One simple alternative
to the functions described in the aforementioned chapter, is to fit a
single polynomial, or a piecewise polynomial (spline) to some given data
points. To avoid a highly fluctuating polynomial, one most often wants
to fit a low-order polynomial to data. This usually means that it is
necessary to fit the polynomial in a least-squares sense, which just is
what the ‘polyfit’ function does.
-- : P = polyfit (X, Y, N)
-- : [P, S] = polyfit (X, Y, N)
-- : [P, S, MU] = polyfit (X, Y, N)
Return the coefficients of a polynomial P(X) of degree N that
minimizes the least-squares-error of the fit to the points ‘[X,
Y]’.
If N is a logical vector, it is used as a mask to selectively force
the corresponding polynomial coefficients to be used or ignored.
The polynomial coefficients are returned in a row vector.
The optional output S is a structure containing the following
fields:
‘R’
Triangular factor R from the QR decomposition.
‘X’
The Vandermonde matrix used to compute the polynomial
coefficients.
‘C’
The unscaled covariance matrix, formally equal to the inverse
of X’*X, but computed in a way minimizing roundoff error
propagation.
‘df’
The degrees of freedom.
‘normr’
The norm of the residuals.
‘yf’
The values of the polynomial for each value of X.
The second output may be used by ‘polyval’ to calculate the
statistical error limits of the predicted values. In particular,
the standard deviation of P coefficients is given by
‘sqrt (diag (s.C)/s.df)*s.normr’.
When the third output, MU, is present the coefficients, P, are
associated with a polynomial in
‘XHAT = (X - MU(1)) / MU(2)’
where MU(1) = mean (X), and MU(2) = std (X).
This linear transformation of X improves the numerical stability of
the fit.
DONTPRINTYET See also: polyval XREFpolyval, *notepolyaffine:
DONTPRINTYET See also: polyval XREFpolyval, polyaffine
XREFpolyaffine, roots XREFroots, vander XREFvander,
zscore XREFzscore.
In situations where a single polynomial isn’t good enough, a solution
is to use several polynomials pieced together. The function ‘splinefit’
fits a piecewise polynomial (spline) to a set of data.
-- : 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)
Fit a piecewise cubic spline with breaks (knots) BREAKS to the
noisy data, X and Y.
X is a vector, and Y is a vector or N-D array. If Y is an N-D
array, then X(j) is matched to Y(:,...,:,j).
P is a positive integer defining the number of intervals along X,
and P+1 is the number of breaks. The number of points in each
interval differ by no more than 1.
The optional property PERIODIC is a logical value which specifies
whether a periodic boundary condition is applied to the spline.
The length of the period is ‘max (BREAKS) - min (BREAKS)’. The
default value is ‘false’.
The optional property ROBUST is a logical value which specifies if
robust fitting is to be applied to reduce the influence of outlying
data points. Three iterations of weighted least squares are
performed. Weights are computed from previous residuals. The
sensitivity of outlier identification is controlled by the property
BETA. The value of BETA is restricted to the range, 0 < BETA < 1.
The default value is BETA = 1/2. Values close to 0 give all data
equal weighting. Increasing values of BETA reduce the influence of
outlying data. Values close to unity may cause instability or rank
deficiency.
The fitted spline is returned as a piecewise polynomial, PP, and
may be evaluated using ‘ppval’.
The splines are constructed of polynomials with degree ORDER. The
default is a cubic, ORDER=3. A spline with P pieces has P+ORDER
degrees of freedom. With periodic boundary conditions the degrees
of freedom are reduced to P.
The optional property, CONSTAINTS, is a structure specifying linear
constraints on the fit. The structure has three fields, "xc",
"yc", and "cc".
"xc"
Vector of the x-locations of the constraints.
"yc"
Constraining values at the locations XC. The default is an
array of zeros.
"cc"
Coefficients (matrix). The default is an array of ones. The
number of rows is limited to the order of the piecewise
polynomials, ORDER.
Constraints are linear combinations of derivatives of order 0 to
ORDER-1 according to
cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).
See also: interp1 XREFinterp1, unmkpp XREFunmkpp,
DONTPRINTYET ppval XREFppval, spline XREFspline, *notepchip:
DONTPRINTYET DONTPRINTYET ppval XREFppval, spline XREFspline, pchip
XREFpchip, ppder XREFppder, ppint XREFppint, *noteDONTPRINTYET DONTPRINTYET ppval XREFppval, spline XREFspline, pchip
XREFpchip, ppder XREFppder, ppint XREFppint,
ppjumps XREFppjumps.
The number of BREAKS (or knots) used to construct the piecewise
polynomial is a significant factor in suppressing the noise present in
the input data, X and Y. This is demonstrated by the example below.
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"})
The piecewise polynomial fit, provided by ‘splinefit’, has continuous
derivatives up to the ORDER-1. For example, a cubic fit has continuous
first and second derivatives. This is demonstrated by the code
## 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"})
When the underlying function to provide a fit to is periodic,
‘splinefit’ is able to apply the boundary conditions needed to manifest
a periodic fit. This is demonstrated by the code below.
## 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"})
More complex constraints may be added as well. For example, the code
below illustrates a periodic fit with values that have been clamped at
the endpoints, and a second periodic fit which is hinged at the
endpoints.
## 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"})
The ‘splinefit’ function also provides the convenience of a ROBUST
fitting, where the effect of outlying data is reduced. In the example
below, three different fits are provided. Two with differing levels of
outlier suppression and a third illustrating the non-robust solution.
## 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
A very specific form of polynomial interpretation is the Padé
approximant. For control systems, a continuous-time delay can be
modeled very simply with the approximant.
-- : [NUM, DEN] = padecoef (T)
-- : [NUM, DEN] = padecoef (T, N)
Compute the Nth-order Padé approximant of the continuous-time delay
T in transfer function form.
The Padé approximant of ‘exp (-sT)’ is defined by the following
equation
Pn(s)
exp (-sT) ~ -------
Qn(s)
Where both Pn(s) and Qn(s) are Nth-order rational functions defined
by the following expressions
N (2N - k)!N! k
Pn(s) = SUM --------------- (-sT)
k=0 (2N)!k!(N - k)!
Qn(s) = Pn(-s)
The inputs T and N must be non-negative numeric scalars. If N is
unspecified it defaults to 1.
The output row vectors NUM and DEN contain the numerator and
denominator coefficients in descending powers of s. Both are
Nth-order polynomials.
For example:
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
The function, ‘ppval’, evaluates the piecewise polynomials, created
by ‘mkpp’ or other means, and ‘unmkpp’ returns detailed information
about the piecewise polynomial.
The following example shows how to combine two linear functions and a
quadratic into one function. Each of these functions is expressed on
adjoined intervals.
x = [-2, -1, 1, 2];
p = [ 0, 1, 0;
1, -2, 1;
0, -1, 1 ];
pp = mkpp (x, p);
xi = linspace (-2, 2, 50);
yi = ppval (pp, xi);
plot (xi, yi);
-- : PP = mkpp (BREAKS, COEFS)
-- : PP = mkpp (BREAKS, COEFS, D)
Construct a piecewise polynomial (pp) structure from sample points
BREAKS and coefficients COEFS.
BREAKS must be a vector of strictly increasing values. The number
of intervals is given by ‘NI = length (BREAKS) - 1’.
When M is the polynomial order COEFS must be of size:
NI-by-(M + 1).
The i-th row of COEFS, ‘COEFS (I,:)’, contains the coefficients for
the polynomial over the I-th interval, ordered from highest (M) to
lowest (0).
COEFS may also be a multi-dimensional array, specifying a
vector-valued or array-valued polynomial. In that case the
polynomial order M is defined by the length of the last dimension
of COEFS. The size of first dimension(s) are given by the scalar
or vector D. If D is not given it is set to ‘1’. In any case
COEFS is reshaped to a 2-D matrix of size ‘[NI*prod(D) M]’.
DONTPRINTYET See also: unmkpp XREFunmkpp, ppval XREFppval, *noteDONTPRINTYET See also: unmkpp XREFunmkpp, ppval XREFppval,
spline XREFspline, pchip XREFpchip, ppder XREFppder,
ppint XREFppint, ppjumps XREFppjumps.
-- : [X, P, N, K, D] = unmkpp (PP)
Extract the components of a piecewise polynomial structure PP.
The components are:
X
Sample points.
P
Polynomial coefficients for points in sample interval. ‘P (I,
:)’ contains the coefficients for the polynomial over interval
I ordered from highest to lowest. If ‘D > 1’, ‘P (R, I, :)’
contains the coefficients for the r-th polynomial defined on
interval I.
N
Number of polynomial pieces.
K
Order of the polynomial plus 1.
D
Number of polynomials defined for each interval.
DONTPRINTYET See also: mkpp XREFmkpp, ppval XREFppval, *noteDONTPRINTYET See also: mkpp XREFmkpp, ppval XREFppval,
spline XREFspline, pchip XREFpchip.
-- : YI = ppval (PP, XI)
Evaluate the piecewise polynomial structure PP at the points XI.
If PP describes a scalar polynomial function, the result is an
array of the same shape as XI. Otherwise, the size of the result
is ‘[pp.dim, length(XI)]’ if XI is a vector, or ‘[pp.dim,
size(XI)]’ if it is a multi-dimensional array.
DONTPRINTYET See also: mkpp XREFmkpp, unmkpp XREFunmkpp, *noteDONTPRINTYET See also: mkpp XREFmkpp, unmkpp XREFunmkpp,
spline XREFspline, pchip XREFpchip.
-- : ppd = ppder (pp)
-- : ppd = ppder (pp, m)
Compute the piecewise M-th derivative of a piecewise polynomial
struct PP.
If M is omitted the first derivative is calculated.
DONTPRINTYET See also: mkpp XREFmkpp, ppval XREFppval, *noteDONTPRINTYET See also: mkpp XREFmkpp, ppval XREFppval,
ppint XREFppint.
-- : PPI = ppint (PP)
-- : PPI = ppint (PP, C)
Compute the integral of the piecewise polynomial struct PP.
C, if given, is the constant of integration.
DONTPRINTYET See also: mkpp XREFmkpp, ppval XREFppval, *noteDONTPRINTYET See also: mkpp XREFmkpp, ppval XREFppval,
ppder XREFppder.
-- : JUMPS = ppjumps (PP)
Evaluate the boundary jumps of a piecewise polynomial.
If there are n intervals, and the dimensionality of PP is d, the
resulting array has dimensions ‘[d, n-1]’.
See also: mkpp XREFmkpp.
Info Catalog
octave: Derivatives / Integrals / Transforms
octave: Polynomial Manipulations
octave: Miscellaneous Functions