d10r9.cpp

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

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

double bessj1(double x)
{
	int tt;
	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 bbb,ccc,aaa,temp,ax,xx,z,y;
    r1 = 72362614232.0;           r2 = -7895059235.0;
    r3 = 242396853.1;            r4 = -2972611.439;
    r5 = 15704.4826;             r6 = -30.16036606;
    s1 = 144725228442.0;          s2 = 2300535178.0;
    s3 = 18583304.74;            s4 = 99447.43394;
    s5 = 376.9991397;            s6 = 1.0;
    p1 = 1.0;                     p2 = 0.00183105;
    p3 = -0.00003516396496;      p4 = 0.000002457520174;
    p5 = -0.000000240337019;
    q1 = 0.04687499995;          q2 = -0.0002002690873;
    q3 = 0.000008449199096;      q4 = -0.00000088228987;
    q5 = 0.000000105787412;
    if (fabs(x) < 8.0 )
	{
        y = x*x;
        aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
        bbb = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
        temp = x * aaa / bbb;
	}
    else
	{
        ax = fabs(x);
        z = 8.0 / ax;
        y = z*z;;
        xx = ax - 2.356194491;
        aaa = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
        bbb = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
        ccc = sqrt(0.636619772 / ax);
		if(x>0) tt=1;
        if(x==0) tt=0;
		if(x>0) tt=-1;
        temp = ccc * (cos(xx) * aaa - z * sin(xx) * bbb * tt);
	}  
	return temp;
}

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  func(double x)
{
 double t;
 t=bessj0( x);
 return t;
}

void zbrak(double x1, double x2, int n, double xb1[], double xb2[], int& nb)
{
	int nbb,i;
	double dx,x,fc,fp;
    nbb = nb;
    nb = 0;
    x = x1;
    dx = (x2 - x1)/ n;
    fp = func(x);
    for (i = 1 ;i<=n;i++)
	{
        x = x + dx;
        fc = func(x);
        if ((fc * fp) < 0.0 )
		{
            nb = nb + 1;
            xb1[nb] = x - dx;
            xb2[nb] = x;
        }
        fp = fc;
        if (nbb == nb)   _c_exit();
    }
}

void  funcd(double x, double &fn,double& df)
{
    fn = bessj0(x);
    df = -bessj1(x);
}

double  rtsafe(double x1, double x2,double xacc)
{
	int maxit,j;
	double temp1,xl,xh,fl,fh,df,swap,temp,dxold,dum,f,dx;
	maxit = 100;
    funcd(x1, fl, df);
    funcd(x2, fh, df);
    if( (fl * fh) >= 0) cout<< "root must be bracketed"<<endl;
    if (fl < 0)
	{
        xl = x1;
        xh = x2;
	}
    else
	{
        xh = x1;
        xl = x2;
        swap = fl;
        fl = fh;
        fh = swap;
	}
    temp1 = 0.5 * (x1 + x2);
    dxold = fabs(x2 - x1);
    dx = dxold;
    funcd(temp1, f, df);
    for (j = 1; j<=maxit; j++)
	{
        dum = (temp1 - xh) * df - f * ((temp1- xl) * df - f);
        if ((dum >= 0)||( fabs(2 * f) > fabs(dxold * df)))
		{
            dxold = dx;
            dx = 0.5 * (xh - xl);
            temp1 = xl + dx;
            if (xl == temp1)  	
			{
				return temp1;
			    _c_exit();
			}
		}
        else
		{
            dxold = dx;
            dx = f / df;
            temp = temp1;
            temp1 = temp1 - dx;
            if (temp == temp1) 	
			{
				return temp1;
				_c_exit();
			}
        }
       if (fabs(dx) < xacc) 	
	   {
		   return temp1;
	      _c_exit();
	   }
       funcd(temp1, f, df);
        if (f < 0 )
		{
            xl = temp1;
            fl = f;
		}
        else
		{
            xh = temp1;
            fh = f;
        }
    }
	return temp1;
    cout<< "rtsafe exceeding maximum iterations."<<endl;
}

void main()
{
    //program d10r9
    //driver for routine rtsafe
	int n,nbmax,nb,i;
	double  x1,x2,xacc,root,xb1[21], xb2[21],ttt;
    n = 100;
    nbmax = 20;
    x1 = 1.0;
    x2 = 50.0;
    nb = nbmax;
    zbrak(x1, x2, n, xb1, xb2, nb);
    cout<<endl;
    cout<<  "Roots of Bessj0:"<<endl;
    cout<<endl;
    cout<< "                    x              f(x)"<<endl;
    for (i = 1; i<=nb; i++)
	{
        xacc = (0.0000010) * (xb1[i] + xb2[i]) / 2.0;
        root =rtsafe(xb1[i], xb2[i], xacc);
        cout<< setw(8)<<"root "<<i;
        cout<< setw(14)<<root;
		ttt=bessj0(root);
        cout<< setw(19)<<ttt<<endl;
    }
}