CubicSpline encapsulates the natural cubic spline of degree 3 with continuity C2. More...

#include <cubicSpline.h>

Collaboration diagram for CubicSpline:

Public Member Functions
	CubicSpline (Int N)
	Constructor.
void	initialize (Real pts[][2])
	Calculates the spline based on the endpoints of the interval of interpolation.
void	dump ()
	Dumps the polynomial coefficients.
Real	operator() (Real v)
	Returns a value approximated by the cubic splice.
Private Attributes
Int	k
	Dimension: k+1 points.
VecReal	x
	The intervals: x[0]..x[1], x[1]..x[2], ...
VecReal	a
VecReal	b
VecReal	c
VecReal	d
	The polynomial coefficients, i = 0,...,k-1.
VecReal	h
VecReal	alp
VecReal	l
VecReal	z
VecReal	mu

Detailed Description

CubicSpline encapsulates the natural cubic spline of degree 3 with continuity C2.

(Natural means that the second derivatives of the spline polynomials are set equal to zero at the endpoints of the interval of interpolation.)

Input: a set of N = k+1 coordinates. Output: a spline as a set of polynomial pieces.

See also:: Algorithm for computing natural cubic splines; Cubic splines module; Mathmatica notebook Natural cubic spline.nb used for tests

Definition at line 25 of file cubicSpline.h.

Constructor & Destructor Documentation

CubicSpline::CubicSpline ( Int N ) [inline]

Constructor.

Only allocates memory for the arrays.

Definition at line 37 of file cubicSpline.h.

        : k( N - 1 )
        , x(N), a(N), b(N), c(N), d(N)
        , h(N), alp(N), l(N), z(N), mu(N)
    {}

Member Function Documentation

void CubicSpline::dump ( ) [inline]

Dumps the polynomial coefficients.

Definition at line 85 of file cubicSpline.h.

References a, b, c, d, k, and x.

Referenced by CubicSpline_sanityCheck().

    {
        std::cout << "Polynomial coefficients:" << std::endl << std::endl;

        for( Int i = 0; i < k; ++i ) {
            std::cout << a[i] << "\t" << b[i] << "\t" << c[i] << "\t" << d[i]
                 << " x-range " << x[i] << " .. " << x[i+1] << std::endl;
        }
    }

void CubicSpline::initialize ( Real pts[][2] ) [inline]

Calculates the spline based on the endpoints of the interval of interpolation.

Definition at line 45 of file cubicSpline.h.

References a, alp, b, c, d, h, k, l, mu, x, and z.

Referenced by CubicSpline_sanityCheck(), and GridUser::cubicSplineShmooth().

    {
        for( Int i = 0; i < k + 1; ++i ) {
            x[i] = pts[i][0];
            a[i] = pts[i][1];
        }

        // Calculate the differences

        for( Int i = 0; i < k; ++i ) {
            h[i] = x[i+1] - x[i];
        }

        for( Int i = 1; i < k; ++i ) {
            alp[i] =  3 * ( a[i+1] - a[i] ) / h[i] - 3 * ( a[i] - a[i-1] ) / h[i-1];
        }

        // Solve the tridiagonal linear equation (Crout's factorization)

        l[0] = 1;  mu[0] = z[0] = 0;

        for( Int i = 1; i < k; ++i ) {
            l[i] = 2 * ( x[i+1] - x[i-1] ) - h[i-1] * mu[i-1];
            mu[i] = h[i] / l[i];
            z[i] = ( alp[i] - h[i-1] * z[i-1] ) / l[i];
        }

        // Compute the polynomial coefficients

        l[k] = 1;  z[k] = c[k] = 0;

        for( Int j = k - 1; j >= 0; --j ) {
            c[j] = z[j] - mu[j] * c[j+1];
            b[j] = ( a[j+1] - a[j] ) / h[j] - h[j] * ( c[j+1] + 2 * c[j] ) / 3;
            d[j] = ( c[j+1] - c[j] ) / ( 3 * h[j] );
        }
    }

Real CubicSpline::operator() ( Real v ) [inline]

Returns a value approximated by the cubic splice.

Definition at line 96 of file cubicSpline.h.

References a, b, c, d, k, and x.

    {
        // Find the polynomial. The polynomial index i goes from 0 to k-1.
        // The interval where the polynomial is valid is x[i-1]..x[i]
        //
        Int i = 0;
        while( i < k - 1 && v >= x[i+1] ) {
            ++i;
        }

        // The interval has not been found. Use the last polynomial for approximation.
        if ( i >= k - 1 ) {
            i = k - 1;
        }

        v -= x[i]; // Subtract the x-value from the reference value

        return a[i] + v * ( b[i] + v * ( c[i] + v * d[i] ) );
    }