三次样条插值算法的C++实现

头文件：

/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                                spline3interp.h
 *
 * Cubic splines interpolation method.
 *
 * For a given points set "Pn" {xn,yn}, this class can find a cubic polynomial
 * in each interval [x_i, x_i+1], such that both of the polynomials and their
 * first order derivative are continuous at the bound of each interval.
 *
 * Zhang Ming, 2010-04, Xi'an Jiaotong University.
 *****************************************************************************/


#ifndef SPLINE3INTERP_H
#define SPLINE3INTERP_H


#include <matrix.h>
#include <interpolation.h>


namespace splab
{

    template <typename Type>
    class Spline3Interp : public Interpolation<Type>
    {

     public:

        using Interpolation<Type>::xi;
        using Interpolation<Type>::yi;

        Spline3Interp( const Vector<Type> &xn, const Vector<Type> &yn,
                       Type d2l=Type(0), Type d2r=Type(0) );
        ~Spline3Interp();

        void calcCoefs();
        Type evaluate( Type x );
        Matrix<Type> getCoefs() const;

    private:

        void derivative2( Vector<Type> &dx, Vector<Type> &d1,
                          Vector<Type> &d2 );

        Type M0,
             Mn;
        Matrix<Type> coefs;

    };
    // class Spline3Interp


    #include <spline3interp-impl.h>

}
// namespace splab


#endif
// SPLINE3INTERP_H

实现文件：

/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                               spline3interp-impl.h
 *
 * Implementation for Spline3Interp class.
 *
 * Zhang Ming, 2010-04, Xi'an Jiaotong University.
 *****************************************************************************/


/**
 * constructors and destructor
 */
template <typename Type>
Spline3Interp<Type>::Spline3Interp( const Vector<Type> &xn,
                                    const Vector<Type> &yn,
                                    Type d2l, Type d2r )
                     : Interpolation<Type>( xn, yn ), M0(d2l), Mn(d2r)
{
}

template <typename Type>
Spline3Interp<Type>::~Spline3Interp()
{
}


/**
 * Compute the second derivative at interpolated points.
 */
template <typename Type>
void Spline3Interp<Type>::derivative2( Vector<Type> &dx, Vector<Type> &d1,
                                  Vector<Type> &d2 )
{
    int N = xi.size(),
        M = N-1;
    Vector<Type> b(M),
                 v(M),
                 y(M),
                 alpha(M),
                 beta(M-1);

	for( int i=1; i<M; ++i )
		b[i] = 2 * (dx[i]+dx[i-1]);

	v[1] = 6*(d1[1]-d1[0]) - dx[0]*d2[0];
	for( int i=1; i<M-1; ++i )
		v[i] = 6 * (d1[i]-d1[i-1]);
	v[M-1] = 6*(d1[M-1]-d1[M-2]) - dx[M-1]*d2[M];

	alpha[1] = b[1];
	for( int i=2; i<M; ++i )
		alpha[i] = b[i] - dx[i]*dx[i-1]/alpha[i-1];

	for( int i=1; i<M-1; ++i )
		beta[i] = dx[i]/alpha[i];

	y[1] = v[1]/alpha[1];
	for( int i=2; i<M; ++i )
		y[i] = (v[i]-dx[i]*y[i-1]) / alpha[i];

	d2[M-1] = y[M-1];
	for( int i=M-2; i>0; --i )
		d2[i] = y[i] - beta[i]*d2[i+1];
}


/**
 * Compute the polynomial' coefsficient in each interval.
 */
template <typename Type>
void Spline3Interp<Type>::calcCoefs()
{
    int N = xi.size(),
        M = N-1;

    Vector<Type> m(N),
                 h(M),
                 d(M);

    m[0] = M0;
    m[M] = Mn;
    for( int i=0; i<M; ++i )
    {
        h[i] = xi[i+1]-xi[i];
        d[i] = (yi[i+1]-yi[i]) / h[i];
    }

    derivative2( h, d, m );

	coefs.resize( M, 4 );
	for( int i=0; i<M; ++i )
	{
		coefs[i][0] = yi[i];
		coefs[i][1] = d[i] - h[i]*(2*m[i]+m[i+1])/6;
		coefs[i][2] = m[i] / 2;
		coefs[i][3] = (m[i+1]-m[i]) / (6*h[i]);
	}
}


/**
 * Compute the value of polynomial at given "x".
 */
template <typename Type>
Type Spline3Interp<Type>::evaluate( Type x )
{
    int k = -1,
        N = xi.size(),
        M = N-1;

	Type dx,
         y;

	for( int i=0; i<M; ++i )
	{
		if( (xi[i]<=x) && (xi[i+1]>=x) )
		{
			k = i;
			dx = x-xi[i];
			break;
		}
	}
	if(k!=-1)
	{
		y = ( ( coefs[k][3]*dx + coefs[k][2] ) * dx + coefs[k][1] ) * dx
		  + coefs[k][0];
		return y;
	}
	else
	{
		cerr << "The value is out of range!" << endl;
		return Type(0);
	}
}


/**
 * Get polynomial' coefsficients.
 */
template <typename Type>
inline Matrix<Type> Spline3Interp<Type>::getCoefs() const
{
	return coefs;
}

测试代码：

/*****************************************************************************
 *                               spline3interp_test.cpp
 *
 * Spline3 interpolation testing.
 *
 * Zhang Ming, 2010-04, Xi'an Jiaotong University.
 *****************************************************************************/


#define BOUNDS_CHECK

#include <iostream>
#include <iomanip>
#include <spline3interp.h>


using namespace std;
using namespace splab;


typedef double   Type;


int main()
{
    // f(x) = 1 / (1+25*x^2)   -1 <= x <= 1
    Type xi[11] = { -1.0,   -0.8,   -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8,   1.0 };
    Type yi[11] = { 0.0385, 0.0588, 0.1,  0.2,  0.5,  1.0, 0.5, 0.2, 0.1, 0.0588, 0.0385 };
    Type Ml = 0.2105, Mr = 0.2105;
    Vector<Type> x( 11, xi ), y( 11, yi );

    Spline3Interp<Type> poly( x, y, Ml, Mr );
    poly.calcCoefs();

    cout << setiosflags(ios::fixed) << setprecision(4);
    cout << "Coefficients of cubic spline interpolated polynomial:   "
         << poly.getCoefs() << endl;

    cout << "The true and interpolated values:" << endl;
    cout << "0.0755" << "   " << poly.evaluate(-0.7) << endl
         << "0.3077" << "   " << poly.evaluate(0.3) << endl
         << "0.0471" << "   " << poly.evaluate(0.9) << endl << endl;

    return 0;
}

运行结果：

Coefficients of cubic spline interpolated polynomial:   size: 10 by 4
0.0385  0.0737  0.1052  0.1694
0.0588  0.1291  0.2069  0.8886
0.1000  0.3185  0.7400  0.8384
0.2000  0.7151  1.2431  13.4078
0.5000  2.8212  9.2877  -54.4695
1.0000  -0.0000 -23.3940        54.4704
0.5000  -2.8212 9.2882  -13.4120
0.2000  -0.7153 1.2410  -0.8225
0.1000  -0.3176 0.7476  -0.9482
0.0588  -0.1323 0.1787  -0.1224

The true and interpolated values:
0.0755   0.0747
0.3077   0.2974
0.0471   0.0472


Process returned 0 (0x0)   execution time : 0.078 s
Press any key to continue.