头文件:
/*
* 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.