d4r24.cpp

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

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>

double bessi1(double x)
{
	double p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9;
	double y,aaa,bbb,ax,temp;
    p1 = 0.5;              p2 = 0.87890594;
    p3 = 0.51498869;       p4 = 0.15084934;
    p5 = 0.02658733;       p6 = 0.00301532;
    p7 = 0.00032411;
    q1 = 0.39894228;       q2 = -0.03988024;
    q3 = -0.00362018;      q4 = 0.00163801;
    q5 = -0.01031555;      q6 = 0.02282967;
    q7 = -0.02895312;      q8 = 0.01787654;
    q9 = -0.00420059;
    if (fabs(x) < 3.75 )
	{
        y = (x / 3.75) *(x / 3.75) ;
        aaa = y * (p4 + y * (p5 + y * (p6 + y * p7)));
        temp = x * (p1 + y * (p2 + y * (p3 + aaa)));
	}
    else
	{
        ax = fabs(x);
        y = 3.75 / ax;
        aaa = exp(ax) / sqrt(ax);
        bbb = y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))));
        temp = aaa * (q1 + y * (q2 + y * (q3 + y * (q4 + bbb))));
	}
    return temp;
}

double bessi0(double x)
{
	double p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9;
	double y,aaa,bbb,ax,temp;
    p1 = 1.0;             p2 = 3.5156229;
    p3 = 3.0899424;       p4 = 1.2067492;
    p5 = 0.2659732;       p6 = 0.0360768;
    p7 = 0.0045813;
    q1 = 0.39894228;      q2 = 0.01328592;
    q3 = 0.00225319;      q4 = -0.00157565;
    q5 = 0.00916281;      q6 = -0.02057706;
    q7 = 0.02635537;      q8 = -0.01647633;
    q9 = 0.00392377;
    if( fabs(x) < 3.75)
	{
        y = (x / 3.75)*(x / 3.75);
        aaa = y * (p5 + y * (p6 + y * p7));
        temp = p1 + y * (p2 + y * (p3 + y * (p4 + aaa)));
	}
    else
	{
        ax = fabs(x);
        y = 3.75 / ax;
        aaa = exp(ax) / sqrt(ax);
        bbb = q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))));
        temp = aaa * (q1 + y * (q2 + y * (q3 + y * bbb)));
	}
	return temp;
}

double bessk1(double x)
{
	double p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7;
	double y,aaa,bbb,ccc,temp; 
    p1 = 1.0;              p2 = 0.15443144;
    p3 = -0.67278579;      p4 = -0.18156897;
    p5 = -0.01919402;      p6 = -0.00110404;
    p7 = -0.00004686;
    q1 = 1.25331414;       q2 = 0.23498619;
    q3 = -0.0365562;       q4 = 0.01504268;
    q5 = -0.00780353;      q6 = 0.00325614;
    q7 = -0.00068245;
    if (x <= 2.0)
	{
        y = x * x / 4;
        aaa = log(x / 2.0) * bessi1(x);
        ccc = y * (p5 + y * (p6 + y * p7));
        bbb = p1 + y * (p2 + y * (p3 + y * (p4 + ccc)));
        temp = aaa + (1.0 / x) * bbb;
	}
    else
	{
        y = (2.0 / x);
        bbb = y * (q5 + y * (q6 + y * q7));
        aaa = q1 + y * (q2 + y * (q3 + y * (q4 + bbb)));
        temp = (exp(-x) / sqrt(x)) * aaa;
	}    
	return temp;
}

double bessk0(double x)
{
	double p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7;
	double y,aaa,bbb,temp;
    p1 = -0.57721566;      p2 = 0.4227842;
    p3 = 0.23069756;       p4 = 0.0348859;
    p5 = 0.00262698;       p6 = 0.0001075;
    p7 = 0.0000074;
    q1 = 1.25331414;       q2 = -0.07832358;
    q3 = 0.02189568;       q4 = -0.01062446;
    q5 = 0.00587872;       q6 = -0.0025154;
    q7 = 0.00053208;
    if( x <= 2.0)
	{
        y = x * x / 4.0;
        bbb = y * (p5 + y * (p6 + y * p7));
        aaa= p1 + y * (p2 + y * (p3 + y * (p4 + bbb)));
        temp = (-log(x / 2.0) * bessi0(x)) + aaa;
	}
    else
	{
        y = (2.0 / x);
        bbb = y * (q5 + y * (q6 + y * q7));
        aaa = q1 + y * (q2 + y * (q3 + y * (q4 + bbb)));
        temp = (exp(-x) / sqrt(x)) * aaa;
	}
	return temp;
}

double bessk(int n,double x)
{
	int j;
	double tox,bkm,bk,bkp,temp;
    if( n < 2)
	{
        cout<<"bad argument n in bessk";
        _c_exit();
    }
    tox = 2.0 / x;
    bkm = bessk0(x);
    bk = bessk1(x);
    for( j = 1;j<=n-1;j++)
	{
        bkp = bkm + j * tox * bk;
        bkm = bk;
        bk = bkp;
    }
    temp = bk;
	return temp;
}

void main()
{
    //program d4r24
    //driver for routine bessk
    int i,n;
	char text[20];
    double nval,value,x; 
    const double pi = 3.1415926;
    fstream fin;
    fin.open("d:\\vc常用数值算法集\\data\\fncval.dat",ios::in);
    while ( strcmp(text,"Kn,")!=0) 
	{
      fin>>text;
	}
	fin>>text;
    fin>>nval; 
    fin>>text;
	cout<<"Modified Bessel Function Kn "<<endl;
    cout<<endl;
    cout<<"     n     x            actual          bessk(n,x)"<<endl;
    for( i = 1; i<=nval; i++)
	{
		fin>>n;
        fin>>x;
		fin>>value;  
		cout<<setprecision(2)<<setw(6)<<n;
		cout<<setprecision(2)<<setw(7)<<x;
        cout<<setprecision(6)<<setw(18)<<value;
		cout<<setprecision(6)<<setw(18)<<bessk(n,x)<<endl;			          
    }
}