nlequations.java

《Java数值计算算法编程》随书代码 已生成ＤＯＣＳ说明文件

	     }
	     
		x.setValue(xx); 
	 }
	
	 /**
	  * 内部函数
	  */
	 private double rnd(Real r)
	 {
	 	int m;
	    double s,u,v,p;
	     
	 	s=65536.0; 
	 	u=2053.0; 
	 	v=13849.0;
	    m=(int)(r.doubleValue()/s); 
	 	r.setValue(r.doubleValue()-m*s);
	    r.setValue(u*r.doubleValue()+v); 
	 	m=(int)(r.doubleValue()/s);
	 	r.setValue(r.doubleValue()-m*s); 
	 	p=r.doubleValue()/s;
	     
	 	return(p);
	 }
	
	 /**
	  * 求实函数或复函数方程一个复根的蒙特卡洛法
	  * 
	  * 调用时，须覆盖计算方程左端函数的模值||f(x, y)||的虚函数
	  *          double func(double x, double y)
	  * 
	  * @param x - 传入初值（猜测解）的实部，返回求得的根的实部
	  * @param y - 传入初值（猜测解）的虚部，返回求得的根的虚部
	  * @param xStart - 均匀分布的端点初值
	  * @param nControlB - 控制参数
	  * @param eps - 控制精度
	  */
	 public void getRootMonteCarlo(Real x, Real y, double xStart, int nControlB, double eps)
	 { 
	     int k;
	     double xx,yy,a,r,z,x1,y1,z1;
	
	 	// 求解条件与初值
	     a=xStart; 
	 	k=1; 
	 	r=1.0; 
	 	xx=x.doubleValue(); 
	 	yy=y.doubleValue();
	     z=func(xx,yy);
	     
	 	// 精度控制求解
	 	while (a>=eps)
	    { 
	     	Real rtmp = new Real(r);
	 		x1=-a+2.0*a*rnd(rtmp); 
	 		r=rtmp.doubleValue();

	 		x1=xx+x1; 
	 		
	     	Real rtmp1 = new Real(r);
	        y1=-a+2.0*a*rnd(rtmp1);
	 		r=rtmp1.doubleValue();
	        
	 		y1=yy+y1;
	         
	 		z1=func(x1,y1);
	         
	 		k=k+1;
	         if (z1>=z)
	         { 
	 			if (k>nControlB) 
	 			{ 
	 				k=1; 
	 				a=a/2.0; 
	 			}
	 		}
	         else
	         { 
	 			k=1; 
	 			xx=x1; 
	 			yy=y1;  
	 			z=z1;
	             if (z<eps)
	             { 
	 				x.setValue(xx); 
	 				y.setValue(yy); 
	 				return; 
	 			}
	         }
	     }
	     
		x.setValue(xx); 
		y.setValue(yy); 
	 }
	
	 /**
	  * 求非线性方程组一组实根的梯度法
	  * 
	  * 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
	  *          double func(double[] x, double[] y)
	  * 
	  * @param n - 方程的个数，也是未知数的个数
	  * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
	  *            返回时存放方程组的一组实根
	  * @param nMaxIt - 迭代次数
	  * @param eps - 控制精度
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootsetGrad(int n, double[] x, int nMaxIt, double eps)
	 { 
	     int l,j;
	     double f,d,s;
	     double[] y = new double[n];
	
	     l=nMaxIt;
	     f=func(x,y);
	
	 	// 控制精度，迭代求解
	     while (f>=eps)
	     { 
	 		l=l-1;
	        if (l==0) 
	 		{ 
	 			return true;
	 		}
	         
	 		d=0.0;
	        for (j=0; j<=n-1; j++) 
	 			d=d+y[j]*y[j];
	        if (d+1.0==1.0) 
	 		{ 
	 			return false;
	 		}
	         
	 		s=f/d;
	        for (j=0; j<=n-1; j++) 
	 			x[j]=x[j]-s*y[j];
	         
	 		f=func(x,y);
	     }
	     
	 	// 是否在有效迭代次数内达到精度
	 	return (nMaxIt>l);
	 }
	
	 /**
	  * 求非线性方程组一组实根的拟牛顿法
	  * 
	  * 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
	  *          double func(double[] x, double[] y)
	  * 
	  * @param n - 方程的个数，也是未知数的个数
	  * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
	  *            返回时存放方程组的一组实根
	  * @param t - 控制h大小的变量，0<t<1
	  * @param h - 增量初值
	  * @param nMaxIt - 迭代次数
	  * @param eps - 控制精度
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootsetNewton(int n, double[] x, double t, double h, int nMaxIt, double eps)
	 { 
	     int i,j,l;
	     double am,z,beta,d;
	
	     double[] y = new double[n];
	 	
	 	 // 构造矩阵
	 	 Matrix mtxCoef = new Matrix(n, n);
	 	 Matrix mtxConst = new Matrix(n, 1);
	     double[] a = mtxCoef.getData();
	     double[] b = mtxConst.getData();
	
	 	 // 迭代求解
	     l=nMaxIt; 
	 	 am=1.0+eps;
	     while (am>=eps)
	     { 
	 		func(x,b);
	         
	 		am=0.0;
	         for (i=0; i<=n-1; i++)
	         { 
	 			z=Math.abs(b[i]);
	             if (z>am) 
	 				am=z;
	         }
	         
	 		if (am>=eps)
	         { 
	 			l=l-1;
	             if (l==0)
	             { 
	                 return false;
	             }
	             
	 			for (j=0; j<=n-1; j++)
	             { 
	 				z=x[j]; 
	 				x[j]=x[j]+h;
	                 
	 				func(x,y);
	                 
	 				for (i=0; i<=n-1; i++) 
	 					a[i*n+j]=y[i];
	                 
	 				x[j]=z;
	             }
	             
	 			// 调用全选主元高斯消元法
	 			LEquations leqs = new LEquations(mtxCoef, mtxConst);
	 			Matrix mtxResult = new Matrix();
	 			if (! leqs.getRootsetGauss(mtxResult))
	 			{
	 				return false;
	 			}
	
	 			mtxConst.setValue(mtxResult);
	 			b = mtxConst.getData();
	
	 			beta=1.0;
	             for (i=0; i<=n-1; i++) 
	 				beta=beta-b[i];
	             
	 			if (beta == 0.0)
	             { 
	                 return false;
	             }
	             
	 			d=h/beta;
	             for (i=0; i<=n-1; i++) 
	 				x[i]=x[i]-d*b[i];
	             
	 			h=t*h;
	         }
	     }
	     
	 	// 是否在有效迭代次数内达到精度
	 	return (nMaxIt>l);
	 }
	
	 /**
	  * 求非线性方程组最小二乘解的广义逆法
	  * 
	  * 调用时，1. 须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
	  *              double func(double[] x, double[] y)
	  *        2. 须覆盖计算雅可比矩阵函数的虚函数
	  *              double FuncMJ(double[] x, double[] y)
	  * 
	  * @param m - 方程的个数
	  * @param n - 未知数的个数
	  * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，要求不全为0，
	  * 			返回时存放方程组的最小二乘解，当m=n时，即是非线性方程组的解
	  * @param eps1 - 最小二乘解的精度控制精度
	  * @param eps2 - 奇异值分解的精度控制精度
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootsetGinv(int m, int n, double[] x, double eps1, double eps2)
	 { 
	 	int i,j,k,l,kk,jt;
        double alpha,z=0,h2,y1,y2,y3,y0,h1;
        double[] p,d,dx;
	     
	 	double[] y = new double[10];
	 	double[] b = new double[10];
	     
	 	// 控制参数
	 	int ka = Math.max(m, n)+1;
	 	double[] w = new double[ka];
	     
	 	// 设定迭代次数为60，迭代求解
	 	l=60; 
	 	alpha=1.0;
	     while (l>0)
	     { 
	 		Matrix mtxP = new Matrix(m, n);
	 		Matrix mtxD = new Matrix(m, 1);
	 		p = mtxP.getData();
	 		d = mtxD.getData();
	
	 		func(x,d);
	        funcMJ(x,p);
	
	 		// 构造线性方程组
	 		LEquations leqs = new LEquations(mtxP, mtxD);
	 		// 临时矩阵
	 		Matrix mtxAP = new Matrix();
	 		Matrix mtxU = new Matrix();
	 		Matrix mtxV = new Matrix();
	 		// 解矩阵
	 		Matrix mtxDX = new Matrix();
	 		// 基于广义逆的最小二乘解
	 		if (! leqs.getRootsetGinv(mtxDX, mtxAP, mtxU, mtxV, eps2))
	 		{ 
	 			return false;
	        }
	         
	 		dx = mtxDX.getData();
	
	 		j=0; 
	 		jt=1; 
	 		h2=0.0;
	        while (jt==1)
	        { 
	 			jt=0;
	             if (j<=2) 
	 				z=alpha+0.01*j;
	             else 
	 				z=h2;
	             
	 			for (i=0; i<=n-1; i++) 
	 				w[i]=x[i]-z*dx[i];
	             
	 			func(w,d);
	             
	 			y1=0.0;
	             for (i=0; i<=m-1; i++) 
	 				y1=y1+d[i]*d[i];
	             for (i=0; i<=n-1; i++)
	 				w[i]=x[i]-(z+0.00001)*dx[i];
	             
	 			func(w,d);
	             
	 			y2=0.0;
	             for (i=0; i<=m-1; i++) 
	 				y2=y2+d[i]*d[i];
	             
	 			y0=(y2-y1)/0.00001;
	             
	 			if (Math.abs(y0)>1.0e-10)
	             { 
	 				h1=y0; h2=z;
	                 if (j==0) 
	 				{ 
	 					y[0]=h1; 
	 					b[0]=h2;
	 				}
	                 else
	                 { 
	 					y[j]=h1; 
	 					kk=0; 
	 					k=0;
	                     while ((kk==0)&&(k<=j-1))
	                     { 
	 						y3=h2-b[k];
	                         if (Math.abs(y3)+1.0==1.0) 
	 							kk=1;
	                         else 
	 							h2=(h1-y[k])/y3;
	                         
	 						k=k+1;
	                     }
	                     
	 					b[j]=h2;
	                     if (kk!=0) 
	 						b[j]=1.0e+35;
	                     
	 					h2=0.0;
	                     for (k=j-1; k>=0; k--)
	 						h2=-y[k]/(b[k+1]+h2);
	                     
	 					h2=h2+b[0];
	                 }
	                 
	 				j=j+1;
	                 if (j<=7) 
	 					jt=1;
	                 else 
	 					z=h2;
	             }
	         }
	         
	 		alpha=z; 
	 		y1=0.0; 
	 		y2=0.0;
	         for (i=0; i<=n-1; i++)
	         { 
	 			dx[i]=-alpha*dx[i];
	             x[i]=x[i]+dx[i];
	             y1=y1+Math.abs(dx[i]);
	             y2=y2+Math.abs(x[i]);
	         }
	         
	 		// 求解成功
	 		if (y1<eps1*y2)
	         { 
	 			return true;
	         }
	         
	 		l=l-1;
	     }
	     
	 	// 求解失败
	 	return false;
	 }
	
	 /**
	  * 求非线性方程组一组实根的蒙特卡洛法
	  * 
	  * 调用时，须覆盖计算方程左端模函数值||F||的虚函数
	  *          double func(int n, double[] x)
	  * 其返回值为Sqr(f1*f1 + f2*f2 + … + fn*fn)
	  * 
	  * @param n - 方程的个数，也是未知数的个数
	  * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
	  *            返回时存放方程组的一组实根
	  * @param xStart - 均匀分布的端点初值
	  * @param nControlB - 控制参数
	  * @param eps - 控制精度
	  */
	 public void getRootsetMonteCarlo(int n, double[] x, double xStart, int nControlB, double eps)
	 { 
	 	 int k,i;
	     double a,r,z,z1;
	
	     double[] y = new double[n];
	     
	 	// 初值
	 	a=xStart; 
	 	k=1; 
	 	r=1.0; 
	
	 	z = func(x);
	
	 	// 用精度控制迭代求解
	     while (a>=eps)
	     { 
	 		for (i=0; i<=n-1; i++)
	 		{
		     	Real rtmp = new Real(r);
	 			y[i]=-a+2.0*a*rnd(rtmp)+x[i];
		 		r=rtmp.doubleValue();
	 		}
	 		z1=func(y);
	         
	 		k=k+1;
	         if (z1>=z)
	         { 
	 			if (k>nControlB) 
	 			{ 
	 				k=1; 
	 				a=a/2.0; 
	 			}
	 		}
	         else
	         { 
	 			k=1; 
	             for (i=0; i<=n-1; i++) 
	 				x[i]=y[i];
	             
	 			// 求解成功
	 			z=z1;
	             if (z<eps)  
	 			{
	 				return;
	 			}
	         }
	     }
	 }
}