| |
- Here we have a code for computing periodic cubic interpolating splines.
These may serve for approximating a given periodic function f
by a two times differentiable function, composed from polynomials of degree not larger than three.
While interpolation by polynomials gives useful results on small intervals only, this type
of approximation is not sensitive against the total length of the interval but only to the distance of the
abscissae.
- The function obtained, we name it S here, is piecewise defined, and
each such piece lives between two successive abscissae of the grid
a = x1 < x2 < ... < xn = b
- These pieces, named Si, i=1,...,n-1 here, are determined from the following conditions:
- Si(x) is a polynomial of degree not greater 3.
- S(x) = Si(x) for xi <= x <= xi+1
- Si(xi) = f(xi) (or yi).
- Si(xi+1) = f(xi+1).
- S'i(xi+1) =
S'i+1(xi+1), i=1,..,n-2 (continuity of S').
- S''i(xi+1) =
S''i+1(xi+1), i=1,...,n-2 (continuity of S'').
Observe that continuity of S itself is obtained by the interpolation condition.
- These conditions lead to a linear system of 4(n-2)+2 equations for
the 4(n-1) coefficients
which is broken down here in a special way into an almost tridiagonal linear system and n-1
uncoupled systems of dimension 3 with a triangular structure.
- This construction is always possible, but not unique.
There remain two degrees of freedom, which allow to fix properties of S
at the boundaries of its interval of definition.
Here we use the periodicity conditions
S'(x1) = S'(xn)
S''(x1) = S''(xn)
which make the solution unique.
- The resultung spline S however will be a periodic function of period
xn-x1 if and only if also
y1 = yn .
We do not check this condition here, hence don't wonder about obscure results if your
input does not satisfy this.
- The resulting linear system has a tridiagonal structure with two additional entries in
the position (1,n) and (n,1) and in the LU decomposition this
makes (only)
the last column of U and the last row of L full. We use a special
Gauss solver for this situation.
- This spline has the remarkable property of minimizing the total curvature integral
∫ab (g''(x))2dx
over the set of all two times continuously differentiable functions g which satisfy the interpolation
and the boundary conditions specified.
- The approximation error of this construction is of the order
h4, with h = max(xi+1-xi),
provided that f is four times continuously differentiable on
[a,b] and periodic with period b-a.
|