d2r1.cpp

一个很实用的采用Neville算法的拉格朗日插值程序。

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

void  polint(double xa[], double ya[], int n, double& x, double& y, double& dy)
{
    double c[11], d[11],dift,dif,ho,hp,w,den;
    int i,m,ns = 1;
    dif = fabs(x - xa[1]);
    for ( i = 1; i<=n; i++)
	{
        dift = fabs(x - xa[i]);
        if (dift < dif)
		{
            ns = i;
            dif = dift;
        }
        c[i] = ya[i];
        d[i] = ya[i];
    }
    y = ya[ns];
    ns = ns - 1;
    for(m = 1; m<=n-1; m++)
	{
        for (i = 1; i<=n-m; i++)
		{
            ho = xa[i] - x;
            hp = xa[i + m] - x;
            w = c[i + 1] - d[i];
            den = ho - hp;
            if (den == 0.0) 
			{
                cout<<"pause"<<endl;
                return;
			}
            den = w / den;
            d[i] = hp * den;
            c[i] = ho * den;
		}
		if (2 * ns < n - m )
            dy = c[ns + 1];
        else
		{
            dy = d[ns];
            ns = ns - 1;
		}
       
        y = y + dy;
	}
}


void main()
{
	//program d2r1
    //river for routine polint
	int nfunc,i;
	int n=9;
	double xa[11], ya[11], x,f,y=0.0,dy=0.0;
    n = 9;    //输入插值节点个数
    const double pi = 3.1415926;
    cout<<"Generation of interpolation tables;"<<endl;
    cout<<" ... sin(x)    0 < x < pi"<<endl;
    cout<<" ... exp(x)    0 < x < 1"<<endl;
    cout<<"How many entries go in these tables?(note: n<10)"<<endl;
    cout<<n<<endl;
    for( nfunc = 1; nfunc<=2; nfunc++)
	{
        if(nfunc == 1 )
		{
            cout<<"sine function from 0 to pi"<<endl;
            for( i = 1; i<=n; i++)
			{
                xa[i] = double(i) * pi / n;
                ya[i] = sin(xa[i]);
            }
		}
        else
		{
			if (nfunc == 2 )
			{
            cout<< "exponential function from 0 to 1"<<endl;
            for ( i = 1; i<=n; i++)
				{
					xa[i] = double(i) * 1 / n;
					ya[i] = exp(xa[i]);
				}
			}
			else
            return;
		}
        cout<<"     x           f(x)   interpolated      error"<<endl;
        for ( i = 1; i<=10; i++)
		{
            if( nfunc == 1 )
			{
                x = (-0.05 + double(i) / 10) * pi;
                f = sin(x);
			}
            else if( nfunc == 2) 
			{
                x = (-0.05 + double(i) / 10);
                f = exp(x);
            }
            polint(xa, ya, n, x, y, dy);
            cout<<setw(10)<<setiosflags(ios::fixed)<<setprecision(6)<<x;
            cout<<setw(12)<<setiosflags(ios::fixed)<<setprecision(6)<<f;
            cout<<setw(12)<<setiosflags(ios::fixed)<<setprecision(6)<<y;
            cout<<setw(16)<<setiosflags(ios::scientific)<<dy<<endl; 
		}
	}    
}