d5r4.cpp

这个是常用的数值算法用VC编写的。相信对大家有用哦。

#include "iostream.h"
#include "math.h"

double func(double x)
{
    return (x * x) * (x * x - 2.0) * sin(x);
}

void chebft(double a, double b, double c[], int n)
{
    int k,nmax = 50;
    double f[51];
    const double pi = 3.14159265358979;
	double bma,bpa,y,sum,fac;
    bma = 0.5 * (b - a);
    bpa = 0.5 * (b + a);
    for (k = 1; k<=n; k++)
	{
        y = cos(pi * (k - 0.5) / n);
        f[k] = func(y * bma + bpa);
    }
    fac = 2.0 / n;
    for (int j = 1; j<=n; j++)
	{
        sum = 0.0;
        for (k = 1; k<=n; k++)
		{
            sum = sum + f[k] * cos((pi * (j - 1)) * ((k - 0.5) / n));
        }
        c[j] = fac * sum;
    }
}

void main()
{
     //program d5r4
     //driver for routine chebft
     int nval = 40;
     double x,y,t0,t1,f,pio2 = 1.5707963;
     double a,b,dum,term,eps = 0.000001;
     double c[41];
     a = -pio2;
     b = pio2;
     chebft(a, b, c, nval);
    //test result
     cout<<endl;
     cout<<"how many terms in chebyshev evaluation?"<<endl;
     //input mval , between 1 and 40, mval=0 to end
     int mval = 20;
     cout<<mval<<endl;
     if ((mval <= 0) || (mval > nval))
	 {
		 return;
	 }
     cout<<"        x            actual     chebyshev fit"<<endl;
	 cout.setf(ios::fixed|ios::right);
	 cout.precision(6);

     for (int i = -8; i<=8; i++)
	 {
         x = i * pio2 / 10.0;
         y = (x - 0.5 * (b + a)) / (0.5 * (b - a));
     //evaluate chebyshev polynomial without using
     //routine chebev
         t0 = 1.0;
         t1 = y;
         f = c[2] * t1 + c[1] * 0.5;
         for (int j = 3; j<=mval; j++)
		 {
             dum = t1;
             t1 = 2.0 * y * t1 - t0;
             t0 = dum;
             term = c[j] * t1;
             f = f + term;
         }
		 cout.width(14);
         cout<<x;
		 cout.width(14);
         cout<<func(x);
		 cout.width(14);
         cout<<f<<endl;
     }
}