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

头文件：

/* * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), [e-mail] * * 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), [e-mail] * * 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 40.0385  0.0737  0.1052  0.16940.0588  0.1291  0.2069  0.88860.1000  0.3185  0.7400  0.83840.2000  0.7151  1.2431  13.40780.5000  2.8212  9.2877  -54.46951.0000  -0.0000 -23.3940        54.47040.5000  -2.8212 9.2882  -13.41200.2000  -0.7153 1.2410  -0.82250.1000  -0.3176 0.7476  -0.94820.0588  -0.1323 0.1787  -0.1224The true and interpolated values:0.0755   0.07470.3077   0.29740.0471   0.0472Process returned 0 (0x0)   execution time : 0.078 sPress any key to continue.