java pnpoly算法_Newton插值算法的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

*/

/*****************************************************************************

* interpolation.h

*

* Base class for interpolation.

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

#ifndef INTERPOLATION_H

#define INTERPOLATION_H

#include

namespace splab

{

template

class Interpolation

{

public:

Interpolation( const Vector &xn, const Vector &yn )

: xi(xn), yi(yn)

{

}

virtual ~Interpolation()

{

}

virtual void calcCoefs() = 0;

virtual Type evaluate( Type x ) = 0;

protected:

Vector xi; // abscissas

Vector yi; // ordinates

};

// class Interpolation

}

// namespace splab

#endif

// INTERPOLATION_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

*/

/*****************************************************************************

* newtoninterp.h

*

* Newton interpolation method.

*

* For a given points set "Pn" {xn,yn}, this class can find a polynomial which

* cross all points of "Pn".

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

#ifndef NEWTONINTERP_H

#define NEWTONINTERP_H

#include

namespace splab

{

template

class NewtonInterp : public Interpolation

{

public:

using Interpolation::xi;

using Interpolation::yi;

NewtonInterp( const Vector &xn, const Vector &yn );

~NewtonInterp();

void calcCoefs();

Type evaluate( Type x );

Vector getCoefs() const;

private:

Vector coefs;

};

// class NewtonInterp

#include

}

// namespace splab

#endif

// NEWTONINTERP_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

*/

/*****************************************************************************

* newtoninterp-impl.h

*

* Implementation for NewtonInterp class.

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

/**

* constructors and destructor

*/

template

NewtonInterp::NewtonInterp( const Vector &xn,

const Vector &yn )

: Interpolation( xn, yn ), coefs(yn)

{

}

template

NewtonInterp::~NewtonInterp()

{

}

/**

* Compute polynomial' coefsficients.

*/

template

void NewtonInterp::calcCoefs()

{

int N = xi.size();

for( int j=1; j

for( int i=N-1; i>=j; --i )

coefs[i] = (coefs[i]-coefs[i-1]) / (xi[i]-xi[i-j]);

}

/**

* Compute the value of polynomial at given "x".

*/

template

Type NewtonInterp::evaluate( Type x )

{

int N = xi.size();

Type y = 0,

tmp = 0;

for( int j=0; j

{

tmp = coefs[j];

for( int i=0; i

tmp *= x-xi[i];

y += tmp;

}

return y;

}

/**

* Get polynomial' coefsficients.

*/

template

inline Vector NewtonInterp::getCoefs() const

{

return coefs;

}

测试代码：

/*****************************************************************************

* newtoninterp_test.cpp

*

* Newton interpolation testing.

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

#define BOUNDS_CHECK

#include

#include

#include

using namespace std;

using namespace splab;

typedef double Type;

int main()

{

Vector x(5),

y(5);

x[0] = 0; x[1] = 30; x[2] = 45; x[3] = 60; x[4] = 90;

y[0] = 0; y[1] = 0.5; y[2] = sqrt(2.0)/2; y[3] = sqrt(3.0)/2; y[4] = 1;

NewtonInterp poly(x,y);

poly.calcCoefs();

cout << setiosflags(ios::fixed) << setprecision(4);

cout << "Coefficients of Newton interpolated polynomial:"

<< poly.getCoefs() << endl;

cout << "The true and interpolated values:" << endl;

cout << sin(10*D2R) << " " << poly.evaluate(10) << endl

<< sin(20*D2R) << " " << poly.evaluate(20) << endl

<< sin(50*D2R) << " " << poly.evaluate(50) << endl << endl;

return 0;

}

运行结果：

Coefficients of Newton interpolated polynomial:size: 5 by 1

0.0000

0.0167

-0.0001

-0.0000

0.0000

The true and interpolated values:

0.1736 0.1734

0.3420 0.3419

0.7660 0.7660

Process returned 0 (0x0) execution time : 0.078 s

Press any key to continue.