d5r6.cpp

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

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

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

double  fint(double x)
{
    //integral of func
	double aaa;
    aaa = 4.0 * x * ((x * x) - 7.0) * sin(x);
    return aaa - ((pow(x , 4) - 14.0 * (x * x) + 28.0)) * cos(x);
}

void chint(double a, double b, double c[], double cint1[], int n)
{
	double con,sum,fac;
    con = 0.25 * (b - a);
    sum = 0.0;
    fac = 1.0;
    for (int j = 2; j<=n - 1; j++)
	{
        cint1[j] = con * (c[j - 1] - c[j + 1]) / (j - 1);
        sum = sum + fac * cint1[j];
        fac = -fac;
    }
    cint1[n] = con * c[n - 1] / (n - 1);
    sum = sum + fac * cint1[n];
    cint1[1] = 2.0 * sum;
}

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;
    }
}

double chebev(double a, double b, double c[], int m, double x)
{
    if ((x - a) * (x - b) > 0.0)
	{
		cout<<"x not in range."<<endl;
	}
    double d = 0.0;
    double sv,y,y2,dd = 0.0;
    y = (2.0 * x - a - b) / (b - a);
    y2 = 2.0 * y;
    for (int j = m; j>=2; j--)
	{
        sv = d;
        d = y2 * d - dd + c[j];
        dd = sv;
    }
    return y * d - dd + 0.5 * c[1];
}

void main()
{
    //program d5r6
    //driver for routine chint
    int nval = 40;
    double a,b,x,aaa,pio2 = 1.5707963;
    double c[41], cint1[41];
    a = -pio2;
    b = pio2;
    chebft(a, b, c, nval);
    //test integral
    cout<<endl;
    cout<<"How many terms in chebyshev evaluation?"<<endl;
    //input mval , between 6 and 40, mval=0 to end
    int mval = 20;
    cout<<mval<<endl;
    if ((mval <= 0) || (mval > nval))
	{
		return;
	}
    chint(a, b, c, cint1, mval);
    cout<<"         x         Actual    Cheby. integ."<<endl;
	cout.setf(ios::fixed|ios::right);
	cout.precision(6);
    for (int i = -8; i<=8; i++)
	{
        x = i * pio2 / 10.0;
		cout.width(13);
        cout<<x;
		cout.width(13);
        cout<<(fint(x) - fint(-pio2));
        aaa = chebev(a, b, cint1, mval, x);
		cout.width(13);
        cout<<aaa<<endl;
    }
}