@c Copyright (C) 1996-2015 John W. Eaton
@c
@c This file is part of Octave.
@c
@c Octave is free software; you can redistribute it and/or modify it
@c under the terms of the GNU General Public License as published by the
@c Free Software Foundation; either version 3 of the License, or (at
@c your option) any later version.
@c
@c Octave is distributed in the hope that it will be useful, but WITHOUT
@c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
@c FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
@c for more details.
@c
@c You should have received a copy of the GNU General Public License
@c along with Octave; see the file COPYING.  If not, see
@c <http://www.gnu.org/licenses/>.

@node Polynomial Manipulations
@chapter Polynomial Manipulations

In Octave, a polynomial is represented by its coefficients (arranged
in descending order).  For example, a vector @var{c} of length
@math{N+1} corresponds to the following polynomial of order
@tex
 $N$
$$
 p (x) = c_1 x^N + \ldots + c_N x + c_{N+1}.
$$
@end tex
@ifnottex
 @var{N}

@example
p(x) = @var{c}(1) x^@var{N} + @dots{} + @var{c}(@var{N}) x + @var{c}(@var{N}+1).
@end example

@end ifnottex

@menu
* Evaluating Polynomials::
* Finding Roots::
* Products of Polynomials::
* Derivatives / Integrals / Transforms::
* Polynomial Interpolation::
* Miscellaneous Functions::
@end menu

@node Evaluating Polynomials
@section Evaluating Polynomials

The value of a polynomial represented by the vector @var{c} can be evaluated
at the point @var{x} very easily, as the following example shows:

@example
@group
N = length (c) - 1;
val = dot (x.^(N:-1:0), c);
@end group
@end example

@noindent
While the above example shows how easy it is to compute the value of a
polynomial, it isn't the most stable algorithm.  With larger polynomials
you should use more elegant algorithms, such as Horner's Method, which
is exactly what the Octave function @code{polyval} does.

In the case where @var{x} is a square matrix, the polynomial given by
@var{c} is still well-defined.  As when @var{x} is a scalar the obvious
implementation is easily expressed in Octave, but also in this case
more elegant algorithms perform better.  The @code{polyvalm} function
provides such an algorithm.

@DOCSTRING(polyval)

@DOCSTRING(polyvalm)

@node Finding Roots
@section Finding Roots

Octave can find the roots of a given polynomial.  This is done by computing
the companion matrix of the polynomial (see the @code{compan} function
for a definition), and then finding its eigenvalues.

@DOCSTRING(roots)

@DOCSTRING(polyeig)

@DOCSTRING(compan)

@DOCSTRING(mpoles)

@node Products of Polynomials
@section Products of Polynomials

@DOCSTRING(conv)

@DOCSTRING(convn)

@DOCSTRING(deconv)

@DOCSTRING(conv2)

@DOCSTRING(polygcd)

@DOCSTRING(residue)

@node Derivatives / Integrals / Transforms
@section Derivatives / Integrals / Transforms

Octave comes with functions for computing the derivative and the integral
of a polynomial.  The functions @code{polyder} and @code{polyint}
both return new polynomials describing the result.  As an example we'll
compute the definite integral of @math{p(x) = x^2 + 1} from 0 to 3.

@example
@group
c = [1, 0, 1];
integral = polyint (c);
area = polyval (integral, 3) - polyval (integral, 0)
@result{} 12
@end group
@end example

@DOCSTRING(polyder)

@DOCSTRING(polyint)

@DOCSTRING(polyaffine)

@node Polynomial Interpolation
@section Polynomial Interpolation

Octave comes with good support for various kinds of interpolation,
most of which are described in @ref{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 @code{polyfit} function does.

@DOCSTRING(polyfit)

In situations where a single polynomial isn't good enough, a solution
is to use several polynomials pieced together.  The function
@code{splinefit} fits a peicewise polynomial (spline) to a set of
data.

@DOCSTRING(splinefit)

The number of @var{breaks} (or knots) used to construct the piecewise
polynomial is a significant factor in suppressing the noise present in
the input data, @var{x} and @var{y}.  This is demonstrated by the example
below.

@example
@group
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"@})
@end group
@end example

@ifnotinfo
@noindent
The result of which can be seen in @ref{fig:splinefit1}.

@float Figure,fig:splinefit1
@center @image{splinefit1,4in}
@caption{Comparison of a fitting a piecewise polynomial with 41 breaks to one
with 11 breaks.  The fit with the large number of breaks exhibits a fast ripple
that is not present in the underlying function.}
@end float
@end ifnotinfo

The piecewise polynomial fit, provided by @code{splinefit}, has
continuous derivatives up to the @var{order}-1.  For example, a cubic fit
has continuous first and second derivatives.  This is demonstrated by
the code

@example
## 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"@})
@end example

@ifnotinfo
@noindent
The result of which can be seen in @ref{fig:splinefit2}.

@float Figure,fig:splinefit2
@center @image{splinefit2,4in}
@caption{Comparison of a piecewise constant, linear, quadratic, cubic, and
quartic polynomials with 8 breaks to noisy data.  The higher order solutions
more accurately represent the underlying function, but come with the
expense of computational complexity.}
@end float
@end ifnotinfo

When the underlying function to provide a fit to is periodic, @code{splinefit}
is able to apply the boundary conditions needed to manifest a periodic fit.
This is demonstrated by the code below.

@example
@group
## 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"@})
@end group
@end example

@ifnotinfo
@noindent
The result of which can be seen in @ref{fig:splinefit3}.

@float Figure,fig:splinefit3
@center @image{splinefit3,4in}
@caption{Comparison of piecewise polynomial fits to a noisy periodic
function with, and without, periodic boundary conditions.}
@end float
@end ifnotinfo

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.

@example
## 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"@})
@end example

@ifnotinfo
@noindent
The result of which can be seen in @ref{fig:splinefit4}.

@float Figure,fig:splinefit4
@center @image{splinefit4,4in}
@caption{Comparison of two periodic piecewise cubic fits to a noisy periodic
signal.  One fit has its endpoints clamped and the second has its endpoints
hinged.}
@end float
@end ifnotinfo

The @code{splinefit} function also provides the convenience of a @var{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.

@example
## 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
@end example

@ifnotinfo
@noindent
The result of which can be seen in @ref{fig:splinefit6}.

@float Figure,fig:splinefit6
@center @image{splinefit6,4in}
@caption{Comparison of two different levels of robust fitting (@var{beta} = 0.25 and 0.75) to noisy data combined with outlying data.  A conventional fit, without
robust fitting (@var{beta} = 0) is also included.}
@end float
@end ifnotinfo

The function, @code{ppval}, evaluates the piecewise polynomials, created
by @code{mkpp} or other means, and @code{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.

@example
@group
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);
@end group
@end example

@DOCSTRING(mkpp)

@DOCSTRING(unmkpp)

@DOCSTRING(ppval)

@DOCSTRING(ppder)

@DOCSTRING(ppint)

@DOCSTRING(ppjumps)

@node Miscellaneous Functions
@section Miscellaneous Functions

@DOCSTRING(poly)

@DOCSTRING(polyout)

@DOCSTRING(polyreduce)