1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
|
#include <iostream>
#include <vector>
using namespace std;
//======================================================================
vector<double> getNewton( const vector<double> &x, vector<double> f ) // Gets Newton polynomial
{
int N = x.size();
vector<double> c(N), temp(N);
c[0] = f[0];
for ( int i = 1; i < N; i++ ) // Compute ith differences
{
for ( int j = 0; j < N - i; j++ ) temp[j] = ( f[j+1] - f[j] ) / ( x[j+i] - x[j] );
f = temp;
c[i] = f[0];
}
return c;
}
//======================================================================
vector<double> standardPolynomial( const vector<double> &c, const vector<double> &x ) // Multiplies out Newton polynomial
{
int N = x.size();
vector<double> a( N, 0.0 ); // a[j] holds coefficient of x^j in final sum
vector<double> p(N), prev(N); // At the ith step, p[ ] is the ith polynomial of Newton factors
p[0] = 1;
a[0] = c[0] * p[0];
for ( int i = 1; i < N; i++ )
{
prev = p;
p[0] = -x[i-1] * prev[0];
a[0] += c[i] * p[0];
for ( int j = 1; j <= i; j++ )
{
p[j] = prev[j-1] - x[i-1] * prev[j];
a[j] += c[i] * p[j];
}
}
return a;
}
//======================================================================
void writeNewtonPolynomial( const vector<double> &c, const vector<double> &x )
{
int N = x.size();
cout << c[0];
for ( int i = 1; i < N; i++ )
{
cout << showpos << c[i];
for ( int j = 0; j < i; j++ ) cout << "(x" << showpos << -x[j] << ")";
}
cout << '\n';
}
//======================================================================
void writeStandardPolynomial( const vector<double> &a )
{
int N = a.size();
cout << a[0];
for ( int i = 1; i < N; i++ ) cout << showpos << a[i] << "x^" << noshowpos << i;
cout << '\n';
}
//======================================================================
double evaluateNewtonPolynomial( const vector<double> &c, const vector<double> &x, double xval )
{
int N = c.size();
double result = c[0], poly = 1.0;
for ( int i = 1; i < N; i++ )
{
poly *= ( xval - x[i-1] );
result += c[i] * poly;
}
return result;
}
//======================================================================
double evaluateStandardPolynomial( const vector<double> &a, double xval )
{
int N = a.size();
double result = a[N-1];
for ( int i = N-2; i >= 0; i-- ) result = a[i] + xval * result;
return result;
}
//======================================================================
int main()
{
vector<double> x = { 1, 2, 3, 4 };
vector<double> y = { 6, 9, 2, 5 };
vector<double> c = getNewton( x, y );
vector<double> a = standardPolynomial( c, x );
cout << "Newton polynomial: "; writeNewtonPolynomial( c, x );
cout << "Standard polynomial: "; writeStandardPolynomial( a );
double test = 2.5;
cout << "Interpolate for x = " << test << '\n';
cout << "Newton: " << evaluateNewtonPolynomial ( c, x, test ) << '\n';
cout << "Standard: " << evaluateStandardPolynomial( a, test ) << '\n';
}
|