d9function.java

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

public class d9Function
{
	double aa,abdevt;
	int ndatat;
	double xt[]=new double[1001];
	double yt[]=new double[1001];
	double arr[]=new double[1001];
	
	void covsrt(double covar[], int ncvm, int ma, int lista[], int mfit)
	{
		int i,j;
		double swap;
		for (j = 1; j<=ma - 1; j++)
		{
			for (i = j + 1; i<=ma; i++)
			{
				covar[(i-1)*ma+j] = 0.0;
			}
	    }
	    for (i = 1; i<=mfit - 1; i++)
		{
	        for (j = i + 1; j<=mfit; j++)
			{
	            if (lista[j] > lista[i])
				{
	                covar[(lista[j]-1)*ma+lista[i]] = covar[(i-1)*ma+j];
				}
	            else
				{
	               covar[(lista[i]-1)*ma+lista[j]] = covar[(i-1)*ma+j];
	            }
	        }
	    }
	    swap = covar[1];;
	    for (j = 1; j<=ma; j++)
		{
	        covar[j] = covar[(j-1)*ma+j];
	        covar[(j-1)*ma+j] = 0.0;
	    }
	    covar[(lista[1]-1)*ma+lista[1]] = swap;
	    for (j = 2; j<=mfit; j++)
		{
	        covar[(lista[j]-1)*ma+lista[j]] = covar[j];
	    }
	    for (j = 2; j<=ma; j++)
		{
	        for (i = 1; i<=j - 1; i++)
			{
	            covar[(i-1)*ma+j] = covar[(j-1)*ma+i];
			}
		}
	}
	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 fit_a;
	double fit_b;
	double fit_siga;
	double fit_sigb;
	double fit_chi2;
	double fit_q;
	void fit(double x[], double y[], int ndata, double sig[], int mwt, double a,
		double b, double siga, double sigb, double chi2, double q)
	{
		int i;
		double sigdat,t,sxoss,wt,ss,sx = 0.0;
		double sy = 0.0;
		double st2 = 0.0;
		b = 0.0;
		if (mwt != 0)
		{
		    ss = 0.0;
		    for (i = 1; i<=ndata; i++)
			{
		        wt = 1.0 / (sig[i] * sig[i]);
		        ss = ss + wt;
		        sx = sx + x[i] * wt;
		        sy = sy + y[i] * wt;
		    }
		}
		else
		{
		    for (i = 1; i<=ndata; i++)
			{
		        sx = sx + x[i];
		        sy = sy + y[i];
		    }
		    ss = ndata;
		}
		sxoss = sx / ss;
		if (mwt != 0)
		{
		    for (i = 1; i<=ndata; i++)
			{
		        t = (x[i] - sxoss) / sig[i];
		        st2 = st2 + t * t;
		        b = b + t * y[i] / sig[i];
		    }
		}
		else
		{
		    for (i = 1; i<=ndata; i++)
			{
		        t = x[i] - sxoss;
		        st2 = st2 + t * t;
		        b = b + t * y[i];
			}
		}
		b = b / st2;
		a = (sy - sx * b) / ss;
		siga = Math.sqrt((1.0 + sx * sx / (ss * st2)) / ss);
		sigb = Math.sqrt(1.0 / st2);
		chi2 = 0.0;
		if (mwt == 0)
		{
		    for (i = 1; i<=ndata; i++)
			{
		        chi2 = chi2 + Math.pow((y[i] - a - b * x[i]) , 2);
		    }
		    q = 1.0;
		    sigdat = Math.sqrt(chi2 / (ndata - 2));
		    siga = siga * sigdat;
		    sigb = sigb * sigdat;
		}
		else
		{
		    for (i = 1; i<=ndata; i++)
			{
		        chi2 = chi2 + Math.pow(((y[i] - a - b * x[i]) / sig[i]) , 2);
			}
		    q = gammq(0.5 * (ndata - 2), 0.5 * chi2);
		}
		fit_a=a;fit_b=b;fit_siga=siga;fit_sigb=sigb;fit_chi2=chi2;fit_q=q;
	}
	void fleg(double x, double pl[], int nl)
	{
		pl[1] = 1.0;
		pl[2] = x;
		if (nl > 2)
		{
			double twox = 2.0 * x;
			double f2 = x;
			double d = 1.0;
			for (int j = 3; j<=nl; j++)
			{
				double f1 = d;
				f2 = f2 + twox;
				d = d + 1.0;
				pl[j] = (f2 * pl[j - 1] - f1 * pl[j - 2]) / d;
			}
		}
	}

	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 lift_chisq;
	void lfit(double x[], double y[], double sig[], int ndata, double a[],int ma,
		  int lista[], int mfit, double covar[], int ncvm, double chisq)
	{
		int i,j,k,kk,ihit;
		double ym,sig2i,wt,sum;
	    double beta[]=new double[51];
	    double afunc[]=new double[51];
	    kk = mfit + 1;
	    for (j = 1; j<=ma; j++)
		{
	        ihit = 0;
	        for (k = 1; k<=mfit; k++)
			{
	            if (lista[k] == j)
				{
					ihit = ihit + 1;
				}
	        }
	        if (ihit == 0)
			{
	            lista[kk] = j;
	            kk = kk + 1;
			}
	        else
			{
				if (ihit > 1)
				{
					System.out.println(" improper set in lista");
					System.exit(1);
				}
	        }
	    }
	    if (kk != (ma + 1))
		{
			System.out.println(" improper set in lista");
			System.exit(1);
		}
	    for (j = 1; j<=mfit; j++)
		{
	        for (k = 1; k<=mfit; k++)
			{
	            covar[(j-1)*ma+k] = 0.0;
	        }
	        beta[j] = 0.0;
	    }
	    for (i = 1; i<=ndata; i++)
		{
	        funcs(x[i], afunc, ma);
	        ym = y[i];
	        if (mfit < ma)
			{
	            for (j = mfit + 1; j<=ma; j++)
				{
	                ym = ym - a[lista[j]] * afunc[lista[j]];
	            }
	        }
	        sig2i = 1.0 / (sig[i] * sig[i]);
	        for (j = 1; j<=mfit; j++)
			{
	            wt = afunc[lista[j]] * sig2i;
	            for (k = 1; k<=j; k++)
				{
	                covar[(j-1)*ma+k] = covar[(j-1)*ma+k] + wt * afunc[lista[k]];
				}
	            beta[j] = beta[j] + ym * wt;
	        }
	    }
	    if (mfit > 1)
		{
	        for (j = 2; j<=mfit; j++)
			{
	            for (k = 1; k<=j - 1; k++)
				{
	                covar[(k-1)*ma+j] = covar[(j-1)*ma+k];
	            }
	        }
	    }
	    gaussj(covar, mfit,ma,beta);
	    for (j = 1; j<=mfit; j++)
		{
	        a[lista[j]] = beta[j];
	    }
	    chisq = 0.0;
	    for (i = 1; i<=ndata; i++)
		{
	        funcs(x[i], afunc, ma);
	        sum = 0.0;
	        for (j = 1; j<=ma; j++)
			{
	            sum = sum + a[j] * afunc[j];
	        }
	        chisq = chisq + ((y[i] - sum) / sig[i]) * ((y[i] - sum) / sig[i]);
	    }
	    covsrt(covar, ncvm, ma, lista, mfit);
		lift_chisq=chisq;
	}

	double medift_a;
	double medift_b;
	double medift_abdev;
	void medfit(double x[], double y[], int ndata, double a, double b, double abdev)
	{
		int j;
	    double bb,del,chisq,f,sx;
	    double sigb,b1,b2,f1,f2,sy = 0.0;
	    double sxy = 0.0;
	    double sxx = 0.0;
		sx=0;
	    for (j = 1; j<=ndata; j++)
		{
	        xt[j] = x[j];
	        yt[j] = y[j];
	        sx = sx + x[j];
	        sy = sy + y[j];
	        sxy = sxy + x[j] * y[j];
	        sxx = sxx + x[j] * x[j];
	    }
	    ndatat = ndata;
	    del = ndata * sxx - sx * sx;
	    aa = (sxx * sy - sx * sxy) / del;
	    bb = (ndata * sxy - sx * sy) / del;
	    chisq = 0.0;
	    for (j = 1; j<=ndata; j++)
		{
	        chisq = chisq + Math.pow((y[j] - (aa + bb * x[j])) , 2);
	    }
	    sigb = Math.sqrt(chisq / del);
	    b1 = bb;
	    f1 = rofunc(b1);
	    b2 = bb + Math.abs(3 * sigb) * sgn(f1);
	    f2 = rofunc(b2);
	    while (f1 * f2 > 0.0)
		{
	        bb = 2.0 * b2 - b1;
	        b1 = b2;
	        f1 = f2;
	        b2 = bb;
	        f2 = rofunc(b2);
		}
	    sigb = 0.01 * sigb;
	    while (Math.abs(b2 - b1) > sigb)
		{
	        bb = 0.5 * (b1 + b2);
	        if (bb == b1 || bb == b2)
			{
				break;
			}
	        f = rofunc(bb);
	        if (f * f1 >= 0.0)
			{
	            f1 = f;
	            b1 = bb;
			}
	        else
			{
	            f2 = f;
	            b2 = bb;
	        }
	    }
	    a = aa;
	    b = bb;
	    abdev = abdevt / ndata;
		medift_a=a;medift_b=b;medift_abdev=abdev;

	}

	
	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 kk,i,j,k,ihit;
	static double ochisq;
	double mrqmin_chisq;
	double mrqmin_alamda;
	static double atry[]=new double[21];
	static double beta[]=new double[21];
	static double da[]=new double[21];
	void mrqmin(double x[], double y[], double sig[], int ndata, double a[],
			int ma, int lista[], int mfit, double covar[], double alpha[],
			int nca, double chisq, double alamda)
	{
		//double ochisq=0;
		if (alamda < 0.0)
		{
		    kk = mfit + 1;
		    for (j = 1; j<=ma; j++)
			{
		        ihit = 0;
		        for (k = 1; k<=mfit; k++)
				{
		            if (lista[k] == j)
					{
						ihit = ihit + 1;
					}
		        }
		        if (ihit == 0)
				{
		            lista[kk] = j;
		            kk = kk + 1;
				}
		        else
				{
					if (ihit > 1)
					{
						System.out.println ("improper permutation in lista");
						System.exit(1);
					}
		        }
		    }
		    if (kk != ma + 1)
			{
				System.out.println ("improper permutation in lista");
				System.exit(1);
			}
		    alamda = 0.001;
		    mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, nca, chisq);
			chisq=mrqcof_chisq;
		    ochisq = chisq;
		    for (j = 1; j<=ma; j++)
			{
		        atry[j] = a[j];
		    }
		}
		for (j = 1; j<=mfit; j++)
		{
		    for (k = 1; k<=mfit; k++)
			{
		        covar[(j-1)*nca+k] = alpha[(j-1)*nca+k];
		    }
		    covar[(j-1)*nca+j] = alpha[(j-1)*nca+j] * (1.0 + alamda);
		    da[j] = beta[j];
		}
		gaussj(covar, mfit, nca, da);
		if (alamda == 0.0)
		{
		    covsrt(covar, nca, ma, lista, mfit);
			for (i=1; i<=20; i++)
			{
				da[i]=0;
				atry[i]=0;
			}
			mrqmin_chisq=chisq;mrqmin_alamda=alamda;
		    return;
		}
		for (j = 1; j<=mfit; j++)
		{
		    atry[lista[j]] = a[lista[j]] + da[j];
		}
		mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, nca, chisq);
		chisq=mrqcof_chisq;
		if (chisq < ochisq)
		{
		    alamda = 0.1 * alamda;
		    ochisq = chisq;
		    for (j = 1; j<=mfit; j++)
			{
		        for (k = 1; k<=mfit; k++)
				{
		            alpha[(j-1)*nca+k] = covar[(j-1)*nca+k];
		        }
		        beta[j] = da[j];
				int temp=lista[j];
		        a[temp] = atry[temp];
		        //a[lista[j]] = atry[lista[j]];
		    }
		}
		else
		{
		    alamda = 10.0 * alamda;
		    chisq = ochisq;
		}
		for (i=1; i<=20; i++)
		{
			da[i]=0;
			atry[i]=0;
		}
		mrqmin_chisq=chisq;	mrqmin_alamda=alamda;
	}

    
   double rofunc(double b)
   {
	    int n1,j,ndata;
		double nml,nmh,d;
		n1 = ndatat + 1;
	    nml = (double)n1 / 2;
	    nmh = n1 - nml;
	    for (j = 1; j<=ndatat; j++)
		{
	        arr[j] = yt[j] - b * xt[j];
	    }
	    sort(ndatat, arr);
	    aa = 0.5 * (arr[(int)nml] + arr[(int)nmh]);
	    double sum = 0.0;
	    abdevt = 0.0;
	    for (j = 1; j<=ndatat; j++)
		{
	        d = yt[j] - (b * xt[j] + aa);
	        abdevt = abdevt + Math.abs(d);
	        sum = sum + xt[j] * sgn(d);
	    }
	    return sum;
	}
	
	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(i=2;i<=ma;i++)
		{
			afunc[i]=x*afunc[i-1];
		}
    }
	
	void gaussj(double A[], int N, int MA,double B[])
	{
		int i,j,k,l,ll,IROW,ICOL;
		double BIG,PIVINV,DUM;
	    int IPIV[]=new int[51];
		int INDXR[]=new int[51];
		int INDXC[]=new int[51];
		ICOL=0;
		IROW=0;
	    for (j=1;j<=N;j++)
		{
	        IPIV[j]=0;
	    }
	    for (i=1;i<=N;i++)
		{
	        BIG=0.0;
	        for (j=1;j<=N;j++)
			{
	            if(IPIV[j]!=1)
				{
	                for(k=1;k<=N;k++)
					{
						if(IPIV[k]==0)
						{
							if(Math.abs(A[(j-1)*MA+k])>=BIG)
							{
								BIG=Math.abs(A[(j-1)*MA+k]);
								IROW=j;
								ICOL=k;
							}
							else if(IPIV[k]>1)
							{
								System.out.println("Singular matrix");
								System.exit(1);
							}
						}
	                }
	            }
	        }
	        IPIV[ICOL]=IPIV[ICOL]+1;
	        if(IROW!=ICOL)
			{
	            for(l=1;l<=N;l++)
				{
	                DUM=(A[(IROW-1)*MA+l]);
	                A[(IROW-1)*MA+l]=A[(ICOL-1)*MA+l];
	                A[(ICOL-1)*MA+l]=DUM;
	            }
	            DUM=B[IROW];
	            B[IROW]=B[ICOL];
	            B[ICOL]=DUM;
			}
	        INDXR[i]=IROW;
	        INDXC[i]=ICOL;
	        if(A[(ICOL-1)*MA+ICOL]==0.0)
			{