d9r4f.java

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

public class d9r4F
{
	void fpoly(double x, double p[], int np)
	{
		p[1] = 1.0;
		for (int j = 2; j <= np; j++)
		{
			p[j] = p[j - 1] * x;
		}
	}
		
	double svdfit_chisq;
	void svdfit(double x[], double y[], double sig[], int ndata, double a[],
    int ma, double u[], double v[], double w[], int mp, int np,
    double chisq, String funcs)
	{
		int i, j;
		double thresh, tmp, tol = 0.00001;
		double b[] = new double[1001];
		double afunc[] = new double[51];
		for (i = 1; i <= ndata; i++)
		{
            if (funcs.compareTo("fpoly") == 0)
			{
				fpoly(x[i], afunc, ma);
			}
            if (funcs.compareTo("fleg") == 0)
			{
				fleg(x[i], afunc, ma);
			}
            tmp = 1.0 / sig[i];
            for (j = 1; j <= ma; j++)
			{
                u[(i - 1) * ma + j] = afunc[j] * tmp;
            }
            b[i] = y[i] * tmp;
		}
		svdcmp(u, ndata, ma, w, v);
		double wmax = 0.0;
		for (j = 1; j <= ma; j++)
		{
            if (w[j] > wmax)
			{
				wmax = w[j];
			}
		}
		thresh = tol * wmax;
		for (j = 1; j <= ma; j++)
		{
            if (w[j] < thresh)
			{
				w[j] = 0.0;
			}
		}
		svbksb(u, w, v, ndata, ma, b, a);
		chisq = 0.0;
		for (i = 1; i <= ndata; i++)
		{
            if (funcs.compareTo("fpoly") == 0)
			{
				fpoly(x[i], afunc, ma);
			}
            if (funcs.compareTo("fleg") == 0)
			{
				fleg(x[i], afunc, ma);
			}
            double sum1 = 0.0;
            for (j = 1; j <= ma; j++)
			{
                sum1 = sum1 + a[j] * afunc[j];
            }
            chisq = chisq + ((y[i] - sum1) / sig[i]) * ((y[i] - sum1) / sig[i]);
		}
		svdfit_chisq = chisq;
	}
  
	void svdvar(double v[], int ma, int np, double w[], double cvm[], int ncvm)
	{
        double wti[] = new double[21];
		int i, j;
		for (i = 1; i <= ma; i++)
		{
			wti[i] = 0.0;
			if (w[i] != 0.0)
			{
				wti[i] = 1.0 / (w[i] * w[i]);
			}
		}
		for (i = 1; i <= ma; i++)
		{
			for (j = 1; j <= i; j++)
			{
				double sum1 = 0.0;
				for (int k = 1; k <= ma; k++)
				{
					sum1 = sum1 + v[(i - 1) * np + k] * v[(j - 1) * np + k] * wti[k];
				}
				cvm[(i - 1) * ma + j] = sum1;
				cvm[(j - 1) * ma + i] = sum1;
			}
		}
	}
	
	void funcs(double x, double afunc[], int ma)
	{
		afunc[1] = 1.0;
		for(int i = 2; i <= ma; i++)
		{
			afunc[i] = x * afunc[i - 1];
		}
    }
	
	
	void svbksb(double u[], double w[], double v[], int m, int n, double b[], double x[])
	{
        double tmp[] = new double[100];
		int i, j, jj;
		double s;
        for (j = 1; j <= n; j++)
		{
            s = 0.0;
            if (w[j] != 0.0)
			{
                for (i = 1; i <= m; i++)
				{
                    s = s + u[(i - 1) * n + j] * b[i];
                }
                s = s / w[j];
            }
            tmp[j] = s;
        }
        for (j = 1; j <= n; j++)
		{
            s = 0.0;
            for (jj = 1; jj <= n; jj++)
			{
                s = s + v[(j - 1) * n + jj] * tmp[jj];
            }
            x[j] = s;
        }
	}

	void svdcmp(double a[],int m,int n,double w[],double v[])
	{
		int nmax=100;
		double rv1[] = new double[nmax];
		int i, j, l, k, nm = 0;
		double g, s, scale, f, h, c, x, y, z; 
		if(m < n)
		{
	 		System.out.println("you must augment a with extra zero rows.");
			System.exit(1);
		}
		g = 0.0;
		scale = 0.0;
		double anorm = 0.0;
		l = 0;
		for (i = 1; i <= n;i++)
		{
	  		l = i + 1;
	  		rv1[i] = scale*g;
	  		g = 0.0;
	  		s = 0.0;
	  		scale = 0.0;
	  		if (i <= m)
			{
                for (k = i; k <= m; k++)
				{
                     scale = scale + Math.abs(a[(k - 1) * n + i]);
				}
                if (scale != 0.0)
				{
                    for (k  =  i; k <= m; k++)
					{
                        a[(k - 1) * n + i] = a[(k - 1) * n + i] / scale;
			       		s = s + a[(k - 1) * n + i] * a[(k - 1) * n + i];
                    }
                    f = a[(i - 1) * n + i];
                    g = -Math.sqrt(s) * sgn(f);
                    h = f * g - s;
                    a[(i - 1) * n + i] = f - g;
                    if (i != n)
					{
                        for (j = l; j <= n; j++)
						{
                            s = 0.0;
				       		for (k = i; k <= m; k++)
							{
                                s = s + a[(k - 1) * n + i] * a[(k - 1) * n + j];
                            }
				       		f = s / h;
	       			   		for (k = i; k <= m; k++)
							{
                                a[(k - 1) * n + j] = a[(k - 1) * n + j] + f * a[(k - 1) * n + i];
		       		   		}
	       	 			}
					}
                    for (k = i; k <= m; k++)
					{
	       				a[(k - 1) * n + i] = scale * a[(k - 1) * n + i];
                    } 
                }
	  		}
	  		w[i] = scale * g;
	  		g = 0.0;
	  		s = 0.0;
	  		scale = 0.0;
	  		if ((i <= m) && (i != n))
			{
                for (k = l; k <= n; k++)
				{
                    scale = scale + Math.abs(a[(i - 1) * n + k]);
                }
                if (scale != 0.0)
				{
                    for (k = l; k <= n; k++)
					{
                        a[(i - 1) * n + k] = a[(i - 1) * n + k] / scale;
						s = s + a[(i - 1) * n + k] * a[(i - 1) * n + k];
                    }
                    f = a[(i - 1) * n + l];
                    g = -Math.sqrt(s) * sgn(f);
                    h = f * g - s;
                    a[(i - 1) * n + l] = f - g;
                    for (k = l; k <= n; k++)
                    {
                        rv1[k] = a[(i - 1) * n + k] / h;
                    }
                    if (i != m)
                    {
                        for (j = l; j <= m; j++)
                        {
                            s = 0.0;
                            for (k = l; k <= n; k++)
                            {
                                s = s + a[(j - 1) * n + k] * a[(i - 1) * n + k];
                            }
                            for (k = l; k <= n; k++)
                            {
                                a[(j - 1) * n + k] = a[(j - 1) * n + k] + s * rv1[k];
                            }
                        }
                    }
                    for (k = l; k <= n; k++)
                    {
                        a[(i - 1) * n + k] = scale * a[(i - 1) * n + k];
                    }
                }
            }
            anorm = Math.max(anorm, (Math.abs(w[i]) + Math.abs(rv1[i])));
        }
        for (i = n; i >= 1; i--)
        {
            if (i < n)
            {     
                if (g != 0.0)
                {
                    for (j = l; j <= n; j++)
                    {
                        v[(j - 1) * n + i] = (a[( i - 1) * n + j] / a[(i - 1) * n + l]) / g;
                    }
                    for (j = l; j <= n; j++)
                    {
                        s = 0.0;
                        for (k = l; k <= n; k++)
                        { 
                            s = s + a[(i - 1) * n + k] * v[(k - 1) * n + j];
                        }
                        for (k = l; k <= n; k++)
                        {
                            v[(k - 1) * n + j] = v[(k - 1) * n + j] + s * v[(k - 1) * n + i];
                        }
                    }
                }
                for (j = l; j <= n; j++)
                {
                    v[(i - 1) * n + j] = 0.0;