d9r8bf.java

解具有导数的多元函数的无约束极小值点或具有导数的非线性方程组的近似解

public class d9r8bF
{
	double fgauss_y; 
	void fgauss(double x, double a[], double y, double dyda[], int na)
	{
		y = 0.0;
        for (int i = 1; i <= na - 1; i = i + 3)
		{
			double arg = (x - a[i + 1]) / a[i + 2];
			double ex = Math.exp(-(arg * arg));
			double fac = a[i] * ex * 2.0 * arg;
			y = y + a[i] * ex;
			dyda[i] = ex;
			dyda[i + 1] = fac / a[i + 2];
			dyda[i + 2] = fac * arg / a[i + 2];
		}
        fgauss_y = y;
	}
		
	double mrqcof_chisq;
	void mrqcof(double x[], double y[], double sig[], int ndata, double a[],
    int ma, int lista[], int mfit, double alpha[], double beta[],
    int nalp, double chisq)
	{
        int i, j, k;
        double wt, ymod, sig2i, dy;
        double  dyda[] = new double[21];
        ymod = 0;
        for (j = 1; j <= mfit; j++)
		{
            for (k = 1; k <= j; k++)
			{
                alpha[(j - 1) * nalp + k] = 0.0;
            }
            beta[j] = 0.0;
        }
        chisq = 0.0;
        for (i = 1; i <= ndata; i++)
		{
            fgauss(x[i], a, ymod, dyda, ma);
            ymod = fgauss_y;
            sig2i = 1.0 / (sig[i] * sig[i]);
            dy = y[i] - ymod;
            for (j = 1; j <= mfit; j++)
			{
                wt = dyda[lista[j]] * sig2i;
                for (k = 1; k <= j; k++)
				{
                    alpha[(j - 1) * nalp + k] = alpha[(j - 1) * nalp + k] + wt * dyda[lista[k]];
                }
                beta[j] = beta[j] + dy * wt;
            }
            chisq = chisq + dy * dy * sig2i;
        }
        for (j = 2; j <= mfit; j++)
		{
            for (k = 1; k <= j - 1; k++)
			{
                alpha[(k - 1) * nalp + j] = alpha[(j - 1) * nalp + k];
            }
        }
        mrqcof_chisq = chisq;
	}

	static int iset;
	static double gset;
	int gasdev_idum;
	double gasdev(int idum)
	{		
        double t, v1, v2, fac, r;
        if (iset == 0.0)
		{
            do
            {
                v1 = 2.0 * ran1(idum) - 1.0;
				idum = ran1_idum;
				v2 = 2.0 * ran1(idum) - 1.0;
				idum = ran1_idum;
                r = v1 * v1 + v2 * v2;
			}while ((r >= 1.0) || (r == 0.0));
			fac = Math.sqrt(-2.0 * Math.log(r) / r);
			gset = v1 * fac;
			t = v2 * fac;
			iset = 1;
			gasdev_idum = idum;
			return t;
		}
        else
		{
			t = gset;
			gasdev_idum = idum;
			iset = 0;
			return t;
        }
	}

	int ran1_idum;
	static double r[] = new double[98];
	static int ix1, ix2, ix3;
	static int iff;
	double ran1(int idum)
	{
        double rm1, rm2, t;
		int m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j;		 
        m1 = 259200; ia1 = 7141; ic1 = 54773; rm1 = 0.0000038580247;
        m2 = 134456; ia2 = 8121; ic2 = 28411; rm2 = 0.0000074373773;
        m3 = 243000; ia3 = 4561; ic3 = 51349;
        if ((idum < 0) || (iff == 0))
		{
            iff = 1;
            ix1 =  (ic1 - idum) % m1;
            ix1 = (ia1 * ix1 + ic1) % m1;
            ix2 = ix1 % m2;
            ix1 = (ia1 * ix1 + ic1) % m1;
            ix3 = ix1 % m3;
            for (j = 1; j <= 97; j++)
			{
                ix1 = (ia1 * ix1 + ic1) % m1;
                ix2 = (ia2 * ix2 + ic2) % m2;
                r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
            }
            idum = 1;
		}
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = (ia2 * ix2 + ic2) % m2;
        ix3 = (ia3 * ix3 + ic3) % m3;
        j = 1 + (int)((97 * ix3) / m3);
        if ((j > 97) || (j < 1)) 
        {
            System.out.println( "abnormal exit");
            System.exit(1);
        }
        t = r[j];
        r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
        ran1_idum = idum;
        return t;
	}
}