d11r8.cpp

共轭梯度数值算法

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

int ncom;
double pcom[51], xicom[51];

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

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 bessj1( double x)
{
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double y,aaa,bbb,ccc,z,ax,xx;
	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))));
		return 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);
		return ccc * (cos(xx) * aaa - z * sin(xx) * bbb * sgn(x));
	}
}

double deriv(double x)
{
    return -bessj1(x);
}


double func2(double x[], int n)
{
    return 1.0 - bessj0(x[1] - 0.5) * bessj0(x[2] - 0.5) * bessj0(x[3] - 0.5);
}


double f1dim(double x)
{
    double xt[51];
    for (int j = 1; j<=ncom; j++)
	{
        xt[j] = pcom[j] + x * xicom[j];
    }
    return func2(xt, ncom);
    //erase xt
}

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

void dfunc(double x[], double df[])
{
    df[1] = bessj1(x[1] - 0.5) * bessj0(x[2] - 0.5) * bessj0(x[3] - 0.5);
    df[2] = bessj0(x[1] - 0.5) * bessj1(x[2] - 0.5) * bessj0(x[3] - 0.5);
    df[3] = bessj0(x[1] - 0.5) * bessj0(x[2] - 0.5) * bessj1(x[3] - 0.5);
}

double dbrent(double ax, double bx, double cx, double tol, double& xmin)
{

    int done,iter,itmax = 100;
    double d,fu,xm,tol1,tol2,a,b,cgold = 0.381966;
    double u,v,w,x,e,fx,fv1,fw,zeps = 0.0000000001;
	double dx,dv,dw,d1,d2,u1,u2,olde,du;
	bool ok1,ok2;
    a = ax;
    if (cx < ax)
	{
		a = cx;
	}
    b = ax;
    if (cx > ax)
	{
		b = cx;
	}
    v = bx;
    w = v;
    x = v;
    e = 0.0;
    fx = func(x);
    fv1 = fx;
    fw = fx;
    dx = deriv(x);
    dv = dx;
    dw = dx;
    for (iter = 1; iter<=itmax; iter++)
	{
        xm = 0.5 * (a + b);
        tol1 = tol * fabs(x) + zeps;
        tol2 = 2.0 * tol1;
		if (fabs(x - xm) <= (tol2 - 0.5 * (b - a)))
		{
			done = -1;
			break;
		}
        else
		{
          done = 0;
		}
        if (fabs(e) > tol1)
		{
            d1 = 2.0 * (b - a);
            d2 = d1;
            if (dw != dx)
			{
				d1 = (w - x) * dx / (dx - dw);
			}
            if (dv != dx)
			{
				d2 = (v - x) * dx / (dx - dv);
			}
            u1 = x + d1;
            u2 = x + d2;
            ok1 = ((a - u1) * (u1 - b) > 0.0) && (dx * d1 <= 0.0);
            ok2 = ((a - u2) * (u2 - b) > 0.0) && (dx * d2 <= 0.0);
            olde = e;
            e = d;
            if (ok1 || ok2)
			{
                if (ok1 && ok2)
				{
                    d = d1;
				}
                else
				{
                    d = d2;
                }
			}
            else
			{
				if (ok1)
				{
					d = d1;
				}
				else
				{
					d = d2;
				}
            }
            if (fabs(d) > fabs(0.5 * olde))
			{
				u = x + d;
				if (u - a < tol2 || b - u < tol2)
				{
					d = fabs(tol1) * sgn(xm - x);
				}
            }
		}
        if (dx >= 0.0)
		{
            e = a - x;
		}
        else
		{
            e = b - x;
        }
        d = 0.5 * e;
        if (fabs(d) >= tol1)
		{
            u = x + d;
            fu = func(u);
		}
        else
		{
            u = x + fabs(tol1) * sgn(d);
            fu = func(u);
            if (fu > fx)
			{
				done = -1;
				break;
			}
            else
			{
				done = 0;
            }
        }
        du = deriv(u);
        if (fu <= fx)
		{
            if (u >= x)
			{
                a = x;
			}
            else
			{
                b = x;
            }
            v = w;
            fv1 = fw;
            dv = dw;
            w = x;
            fw = fx;
            dw = dx;
            x = u;
            fx = fu;
            dx = du;
		}
        else
		{
            if (u < x)
			{
                a = u;
			}
            else
			{
                b = u;
			}
            if (fu <= fw || w == x)
			{
                v = w;
                fv1 = fw;
                dv = dw;
                w = u;
                fw = fu;
                dw = du;
			}
            else
			{
				if (fu <= fv1 || v == x || v == w)
				{
					v = u;
					fv1 = fu;
					dv = du;
				}
			}
		}
	}	
    if (!done)
	{
        cout<<"dbrent exceeded maximum iterations."<<endl;
		exit(1);
	}
    else
	{
        xmin = x;
        return fx;
    }
}

void linmin(double p[], double xi[], int n, double& fret)
{
	int j;
    double tol = 0.0001;
    ncom = n;
    for (j = 1; j<=n; j++)
	{
        pcom[j] = p[j];
        xicom[j] = xi[j];
    }
    double fa,fx,fb,bx,ax = 0.0;
    double xmin,xx = 1.0;
    mnbrak(ax, xx, bx, fa, fx, fb);
    fret = dbrent(ax, xx, bx, tol, xmin);
    for (j = 1; j<=n; j++)
	{
        xi[j] = xmin * xi[j];
        p[j] = p[j] + xi[j];
    }
}


void frprmn(double p[], int n, double ftol, int& iter, double& fret)
{
    int j,its,itmax = 200;
    double eps = 0.0000000001;
    double g[551], h[51], xi[51];
    double fp = func2(p, n);
    dfunc(p, xi);
    for (j = 1; j<=n; j++)
	{
        g[j] = -xi[j];
        h[j] = g[j];
        xi[j] = h[j];
    }
    for (its = 1; its<=itmax; its++)
	{
        iter = its;
        linmin(p, xi, n, fret);
        if (2.0 * fabs(fret - fp) <= ftol * (fabs(fret) + fabs(fp) + eps))
		{
            break;
        }
        fp = func2(p, n);
        dfunc(p, xi);
        double gg = 0.0;
        double dgg = 0.0;
        for (j = 1; j<=n; j++)
		{
            gg = gg + g[j] * g[j];
            //dgg = dgg + xi(j) ^ 2              polak-ribiere 法
            dgg = dgg + (xi[j] + g[j]) * xi[j];  //fletcher-reeves 法
        }
        if (gg == 0.0)
		{
			break;
		}
        double gam = dgg / gg;
        for (j = 1; j<=n; j++)
		{
            g[j] = -xi[j];
            h[j] = g[j] + gam * h[j];
            xi[j] = h[j];
        }
    }
    if (iter > itmax)
	{
		cout<<"  frpr maximum iterations exceeded";
		cout<<endl;
	}
}

void main()
{
    //program d11r8
    //driver for routine frprmn
    int iter,k,ndim = 3;
    double fret,angl,pio2 = 1.5707963;
    double ftol = 0.000001;
    double p[4];
    cout<<"Program finds the minimum of a function"<<endl;
    cout<<"with different trial starting vectors."<<endl;
    cout<<"True minimum is (0.5, 0.5, 0.5)"<<endl;
	cout.setf(ios::fixed|ios::left);
	cout.precision(6);
    for (k = 0; k<=4; k++)
	{
        angl = pio2 * k / 4.0;
        p[1] = 2.0 * cos(angl);
        p[2] = 2.0 * sin(angl);
        p[3] = 0.0;
        cout<<endl;
        cout<<"Starting vector: (";
		cout.width(12);
		cout<<p[1];
		cout.width(12);
        cout<<p[2];
		cout.width(7);
		cout<<p[3]<<")"<<endl;
        frprmn(p, ndim, ftol, iter, fret);
        cout<<"Iterations: "<<iter<<endl;
        cout<<"Solution vector: (";
		cout.width(12);
		cout<<p[1];
		cout.width(12);
        cout<<p[2];
		cout.width(7);
		cout<<p[3]<<")"<<endl;
		cout.setf(ios::scientific);
        cout<<"func. value at solution "<<fret<<endl;
		cout.unsetf(ios::scientific);
		cout.setf(ios::fixed);
		cout.precision(6);
    }
}