d10r2.cpp

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

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

double bessj0(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double y,bbb,ccc,aaa,temp,eee,ddd,ax,xx,z;
    p1 = 1.0;                  p2 = -0.001098628627;
    p3 = 0.00002734510407;    p4 = -0.000002073370639;
    p5 = 2.093887211e-07;
    q1 = -0.01562499995;      q2 = 0.0001430488765;
    q3 = -0.000006911147651;  q4 = 7.621095161e-07;
    q5 = -9.34945152e-08;
    r1 = 57568490574.0;        r2 = -13362590354.0;
    r3 = 651619640.7;         r4 = -11214424.18;
    r5 = 77392.33017;         r6 = -184.9052456;
    s1 = 57568490411.0;        s2 = 1029532985.0;
    s3 = 9494680.718;         s4 = 59272.64853;
    s5 = 267.8532712;        s6 = 1.0;
    if (fabs(x) < 8.0)
	{
       y = x * x;
       bbb = y * (r4 + y * (r5 + y * r6));
       aaa = r1 + y * (r2 + y * (r3 + bbb));
       ccc = y * (s3 + y * (s4 + y * (s5 + y * s6)));
       temp = aaa / (s1 + y * (s2 + ccc));
	}
    else
	{
       ax = fabs(x);
       z = 8.0 / ax;
       y = z * z;
       xx = ax - 0.785398164;
       ccc = y * (p3 + y * (p4 + y * p5));
       aaa = p1 + y * (p2 + ccc);
       ddd = y * (q3 + y * (q4 + y * q5));
       eee = z * sin(xx) * (q1 + y * (q2 + ddd));
       temp = sqrt(0.636619772 / ax) * (cos(xx) * aaa - eee);
	}
	return temp;
}

double fun(double x)
{
 double t;
 t=bessj0( x);
 return t;
}

void zbrac(double &x1, double &x2, bool &succes)
{
	int j,ntry;
	double factor,f1,f2;
    factor = 1.6;
    ntry = 50;
    if (x1 == x2) cout<< "you have to guess an initial range"<<endl;
    f1 = fun(x1);
    f2 = fun(x2);
    succes = true;
    for (j = 1; j<=ntry; j++)
	{
        if ((f1 * f2) < 0.0) 
		goto loop;
        if (fabs(f1) < fabs(f2))
		{
            x1 = x1 + factor * (x1 - x2);
            f1 = fun(x1);
		}
        else
		{
            x2 = x2 + factor * (x2 - x1);
            f2 = fun(x2);
        }
	}
    succes = false;
loop:cout<<endl;
}

void main()
{
	//program d10r2
    //driver for routine zbrac
	bool succes;
	int i;
	double x1,x2;
    cout<<endl;
    cout<<"  bracketing values:      function values:"<<endl;
    cout<<endl;
    cout<<"     x1      x2      bessj0(x1)    bessj0(x2)";
	cout<<setprecision(2)<<setiosflags(ios::fixed);
    for (i = 1; i<=10; i++)
	{
        x1 = double(i);
        x2 = x1 + 1;
        zbrac(x1, x2, succes);
        if (succes)
		{
            cout<<setprecision(2)<<setw(8)<<x1;
            cout<<setprecision(2)<<setw(8)<<x2;
            cout<<setprecision(6)<<setw(14)<<bessj0(x1);
            cout<<setprecision(6)<<setw(14)<<bessj0(x2);
        }
    }
	cout<<endl;
}