d2r8.cpp

数值计算c++源代码,包括各种算法。很有用的。

# 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 polcof(double xa[],double ya[],int n, double cof[])
{
	double   x[15], y[15],xmin,dy,kk;
	int j,i,k;
	kk=0.0;
    for (j = 1;j<=n;j++)
	{
        x[j] = xa[j];
        y[j] = ya[j];
    }
    for (j = 1;j<=n;j++)
	{
        polint(x, y, n + 1 - j, kk, cof[j], dy);
        xmin = 1e+38;
        k = 0;
        for( i = 1; i<=n + 1 - j; i++)
		{
			
            if (fabs(x[i]) < xmin )
			{
                xmin = fabs(x[i]);
                k = i;
            }
            if (x[i] !=0.0)
				y[i] = (y[i] - cof[j]) / x[i];
        }
        if (k < n + 1 - j)
		{
            for (i = k + 1; i<=n + 1 - j; i++)
			{
                y[i - 1] = y[i];
                x[i - 1] = x[i];
            }
        }
    }
}
   
void main()
{
    //program d2r8
    //driver for routine polcof
    int nfunc,i,j,np = 5;
    double x,f,sum,xa[6], ya[6], coeff[6],pi = 3.141593;
    for (nfunc = 1; nfunc<=2; nfunc++)
	{
        if (nfunc == 1)
		{
            cout<<endl;
            cout<<"sine function from 0 to pi"<<endl;
            for (i = 1; i<=np; i++)
			{
                xa[i] = i * pi /double(np);
                ya[i] = sin(xa[i]);
            }
		}
        else 
		{
			if (nfunc ==2)
			{
            cout<<endl;
            cout<< "Exponential function from 0 to 1"<<endl;
            for (i = 1;i<=np;i++)
				{
					xa[i] = 1.0 * i / double(np);
					ya[i] = exp(xa[i]);
				}
			}
            else
            return;
        }
        polcof(xa, ya, np, coeff);
        cout<<"coefficients"<<endl;
        for (i = 1;i<=np;i++)
            cout<<setw(12)<<coeff[i];
		cout<<endl;
        cout<< "           x          f[x]       polynomial"<<endl;
        for (i = 1; i<=10; i++)
		{
            if (nfunc == 1 )
			{
                x = (-0.05 + i / 10.0) * pi;
                f = sin(x);
			}
            else 
			{
				if (nfunc == 2)
				{
					x = (-0.05 + double(i) / 10.0);
					f = exp(x);
				}
            }
            sum = coeff[np];
            for (j = np - 1; j>=1; j--)
				{
					sum = coeff[j] + sum * x;
				}
            cout<< setw(14)<<x;
            cout<< setw(14)<<f;
            cout<< setw(14)<<sum<<endl;
		}
        cout<<  "**********************************"<<endl;
    }
}