d11r5.cpp

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

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

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

double famoeb(double x[])
{
	double aaa;
    aaa = (x[1] - 0.5) * (x[1] - 0.5) + (x[2] - 0.6) * (x[2] - 0.6) + (x[3] - 0.7) * (x[3] - 0.7);
    return 0.6 - bessj0(aaa);
}

void erase(double pbar[], double prr[], double pr[])
{
	for (int i=0; i<=20; i++)
	{
		pbar[i]=0;
		prr[i]=0;
		pr[i]=0;
	}
}

void amoeba(double p[], double y[], int mp, int np, int ndim, double ftol, int& iter)
{
    int i,j,ihi,inhi,mpts,nmax = 20;
    double ypr,yprr,rtol,alpha = 1.0;
    double beta = 0.5;
    double gamma = 2.0;
    int itmax = 500;
    double pr[21], prr[21], pbar[21];
    mpts = ndim + 1;
    iter = 0;
    do
	{
		int ilo = 1;
        if (y[1] > y[2])
		{
            ihi = 1;
            inhi = 2;
		}
        else
		{
            ihi = 2;
            inhi = 1;
		}
        for (i = 1; i<=mpts; i++)
		{
            if (y[i] < y[ilo])
			{
			    ilo = i;
			}
            if (y[i] > y[ihi])
			{
                inhi = ihi;
                ihi = i;
			}
            else
			{
			    if (y[i] > y[inhi])
				{
				    if (i != ihi)
					{
					    inhi = i;
					}
				}
			}
		}
        rtol = 2.0 * fabs(y[ihi] - y[ilo]) / (fabs(y[ihi]) + fabs(y[ilo]));
        if (rtol < ftol)
		{
		    erase(pbar, prr, pr);
		    return;
		}
        if (iter == itmax)
		{
            cout<<" amoeba exceeding maximum iterations."<<endl;
            return;
		}
        iter = iter + 1;
        for (j = 1; j<=ndim; j++)
		{
            pbar[j] = 0.0;
		}
        for (i = 1; i<=mpts; i++)
		{
            if (i != ihi)
			{
                for (j = 1; j<=ndim; j++)
				{
                    pbar[j] = pbar[j] + p[(i-1)*np+j];
				}
			}
		}
        for (j = 1; j<=ndim; j++)
		{
            pbar[j] = pbar[j] / ndim;
            pr[j] = (1.0 + alpha) * pbar[j] - alpha * p[(ihi-1)*np+j];
		}
        ypr = famoeb(pr);
        if (ypr <= y[ilo])
		{
            for (j = 1; j<=ndim; j++)
			{
                prr[j] = gamma * pr[j] + (1.0 - gamma) * pbar[j];
			}
            yprr = famoeb(prr);
            if (yprr < y[ilo])
			{
                for (j = 1; j<=ndim; j++)
				{
                    p[(ihi-1)*np+j] = prr[j];
				}
                y[ihi] = yprr;
			}
            else
			{
                for (j = 1; j<=ndim; j++)
				{
                    p[(ihi-1)*np+j] = pr[j];
				}
                y[ihi] = ypr;
			}
		}
        else
		{
		    if (ypr >= y[inhi])
			{
			    if (ypr < y[ihi])
				{
				    for (j = 1; j<=ndim; j++)
					{
						p[(ihi-1)*np+j] = pr[j];
					}
				}
				y[ihi] = ypr;
			
				for (j = 1; j<=ndim; j++)
				{
					prr[j] = beta * p[(ihi-1)*np+j] + (1.0 - beta) * pbar[j];
				}
				yprr = famoeb(prr);
				if (yprr < y[ihi])
				{
					for (j = 1; j<=ndim; j++)
					{
						p[(ihi-1)*np+j] = prr[j];
					}
					y[ihi] = yprr;
				}

				else
				{
					for (i = 1; i<=mpts; i++)
					{
						if (i != ilo)
						{
							for (j = 1; j<=ndim; j++)
							{
								pr[j] = 0.5 * (p[(i-1)*np+j] + p[(ilo-1)*np+j]);
								p[(i-1)*np+j] = pr[j];
							}
							y[i] = famoeb(pr);
						}
					}
				}
			}
			else
			{
				for (j = 1; j<=ndim; j++)
				{
					p[(ihi-1)*np+j] = pr[j];
				}
				y[ihi] = ypr;
			}
		}
    }while(1);
}

void main()
{
    //program d11r5
    //driver for routine amoeba
    int iter,ndim,np = 3;
    int i,j,mp = 4;
    double ftol = 0.000001;
    double p[13], x[4], y[5];
    p[1] = 0.0;  p[2] = 0.0;  p[3] = 0.0;
    p[4] = 1.0;  p[5] = 0.0;  p[6] = 0.0;
    p[7] = 0.0;  p[8] = 1.0;  p[9] = 0.0;
    p[10] = 0.0; p[11] = 0.0; p[12] = 1.0;
    ndim = np;
    for (i = 1; i<=mp; i++)
	{
        for (j = 1; j<=np; j++)
		{
            x[j] = p[(i-1)*np+j];
        }
        y[i] = famoeb(x);
    }
    amoeba(p, y, mp, np, ndim, ftol, iter);
    cout<<endl;
    cout<<"Iterations: "<<iter<<endl;
    cout<<"Vertices of final 3-D simplex and"<<endl;
    cout<<"function values at the vertices:"<<endl;
    cout<<"i       x(i)         y(i)       z(i)       function"<<endl;
	cout.setf(ios::fixed|ios::left);
	cout.precision(6);
    for (i = 1; i<=mp; i++)
	{
		cout.width(6);
        cout<<i;
        for (j = 1; j<=np; j++)
		{
			cout.width(12);
			cout<<p[(i-1)*np+j];
        }
		cout.width(12);
        cout<<y[i]<<endl;
    }
	cout<<endl;
    cout<<"True minimum is at ( 0.5, 0.6, 0.7)"<<endl;
}