d11r2.cpp

使用VC++编写的大量数学算法的源代码

#include "iostream.h"
#include "math.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;
}

int sgn(double x)
{
	int temp;
	if(x>0) temp=1;
	if(x=0) temp=0;
	if(x<0) temp=-1;
	return temp;
}

double bessj1(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 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);
        temp = ccc * (cos(xx) * aaa - z * sin(xx) * bbb * sgn(x));
	}  
	return temp;
}

double func(double x)
{
    return bessj0(x);
}

void mnbrak(double& ax, double& bx, double& cx, double& fa, double& fb, double& fc)
{
    double r,q,dum,gold = 1.618034;
    int glimit = 100;
    double u,ulim,fu,tiny = 1e-20;
    fa = func(ax);
    fb = func(bx);
    if (fb > fa)
	{
        dum = ax;
        ax = bx;
        bx = dum;
        dum = fb;
        fb = fa;
        fa = dum;
    }
    cx = bx + gold * (bx - ax);
    fc = func(cx);
	while (fb >= fc)
	{
        r = (bx - ax) * (fb - fc);
        q = (bx - cx) * (fb - fa);
        dum = q - r;
        if (fabs(dum) < tiny)
		{
			dum = tiny;
		}
        u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
        ulim = bx + glimit * (cx - bx);
        if ((bx - u) * (u - cx) > 0)
		{
            fu = func(u);
            if (fu < fc)
			{
                ax = bx;
                fa = fb;
                bx = u;
                fb = fu;
                return;
			}
            else
			{
				if (fu > fb)
				{
					cx = u;
					fc = fu;
					return;
				}
            }
            u = cx + gold * (cx - bx);
            fu = func(u);
		}
        else
		{
			if ((cx - u) * (u - ulim) > 0)
			{
				fu = func(u);
				if (fu < fc)
				{
					bx = cx;
					cx = u;
					u = cx + gold * (cx - bx);
					fb = fc;
					fc = fu;
					fu = func(u);
				}
            }
			else
			{
				if ((u - ulim) * (ulim - cx) >= 0)
				{
					u = ulim;
					fu = func(u);
				}
				else
				{
					u = cx + gold * (cx - bx);
					fu = func(u);
				}
			}
		}
		ax = bx;
		bx = cx;
		cx = u;
		fa = fb;
		fb = fc;
		fc = fu;
	}
}

double golden(double ax, double bx, double cx, double tol, double& xmin)
{
    double x0,x1,x2,x3,f0,f1,f2,f3,r= 0.61803399;
    double c = 0.38196601;
    x0 = ax;
    x3 = cx;
    if (fabs(cx - bx) > fabs(bx - ax))
	{
		x1 = bx;
        x2 = bx + c * (cx - bx);
	}
    else
	{
        x2 = bx;
        x1 = bx - c * (bx - ax);
    }
    f1 = func(x1);
    f2 = func(x2);
    while (fabs(x3 - x0) > tol * (fabs(x1) + fabs(x2)))
	{
        if (f2 < f1)
		{
            x0 = x1;
            x1 = x2;
            x2 = r * x1 + c * x3;
            f0 = f1;
            f1 = f2;
            f2 = func(x2);
		}
        else
		{
            x3 = x2;
            x2 = x1;
            x1 = r * x2 + c * x0;
            f3 = f2;
            f2 = f1;
            f1 = func(x1);
        }
    }
    if (f1 < f2)
	{
        xmin = x1;
        return f1;
	}
    else
	{
        xmin = x2;
        return f2;
    }
}

void main()
{
     //program d11r2
     //driver for routine golden
     double g,ax,bx,cx,fa,fb,fc,tol = 0.000001;
     double xmin,eql = 0.001;
     double amin[21];
     int iflag,i,nmin = 0;
     cout<<endl;
     cout<<"minima of the function Bessj0"<<endl;
     cout<<endl;
     cout<<"  min. #     x         Bessj0(x)     Bessj1(x)"<<endl;
	 //cout.setf(ios::fixed|ios::left);
	 //cout.precision(6);
     for (i = 1; i<=100; i++)
	 {
         ax = i;
         bx = i + 1.0;
         mnbrak(ax, bx, cx, fa, fb, fc);
         g = golden(ax, bx, cx, tol, xmin);
         if (nmin == 0)
		 {
             amin[1] = xmin;
             nmin = 1;
			 cout.width(4);			 cout<<nmin;
			 cout.width(13);		 cout<<xmin;
			 cout.width(15);		 cout<<bessj0(xmin);
			 cout.width(16);		 cout<<bessj1(xmin);
			 cout<<endl;
		 }
         else
		 {
             iflag = 0;
             for (int j = 1; j<=nmin; j++)
			 {
                 if (fabs(xmin - amin[j]) <= eql * xmin)
				 {
					 iflag = 1;
				 }
             }
             if (iflag == 0)
			 {
                 nmin = nmin + 1;
                 amin[nmin] = xmin;
				 cout.width(4);		 cout<<nmin;
				 cout.width(13);	 cout<<xmin;
				 cout.width(15);	 cout<<bessj0(xmin);
				 cout.width(16);	 cout<<bessj1(xmin);
				 cout<<endl;
             }
         }
     }
}