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.

    {
        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.

    {
        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.

    {
        // 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] ) );
    }