nlequations.java

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

/*
 * 求解非线性方程组的类 NLEquations
 * 
 * 周长发编制
 */
package javaalgorithm.algorithm;

/**
 * 求解非线性方程组的类 NLEquations
 *
 * @author 周长发
 * @version 1.0
 */
public abstract class NLEquations 
{
	/**
	 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
	 * 
	 * @param x - 变量
	 * @return 函数值
	 */
	protected double func(double x)
	{
		return 0;
	}
	
	/**
	 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
	 * 
	 * @param x - 变量值数组
	 * @return 函数值
	 */
	protected double func(double[] x)
	{
		return 0;
	}


	/**
	 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
	 * 
	 * @param x - 变量
	 * @param y - 函数值数组
	 */
	protected void func(double x, double[] y)
	{
	}
	
	/**
	 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
	 * 
	 * @param x - 二元函数的变量
	 * @param y - 二元函数的变量
	 * @return 函数值
	 */
	protected double func(double x, double y)
	{
		return 0;
	}
	
	/**
	 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
	 * 
	 * @param x - 二元函数的变量值数组
	 * @param y - 二元函数的变量值数组
	 * @return 函数值
	 */
	protected double func(double[] x, double[] y)
	{
		return 0;
	}
	
	/**
	 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
	 * 
	 * @param x - 已知变量值数组
	 * @param p - 已知函数值数组
	 */
	protected void funcMJ(double[] x, double[] p)
	{
	}

	/**
     * 基本构造函数
     */
	 public NLEquations()
	 {
	 }

	 /**
	  * 求非线性方程实根的对分法
	  * 
	  * 调用时，须覆盖计算方程左端函数f(x)值的虚函数
	  *        double func(double x)
	  * 
	  * @param nNumRoots - 在[xStart, xEnd]内实根个数的预估值
	  * @param x - 一维数组，长度为m。返回在区间[xStart, xEnd]内搜索到的实根，
	  *            实根个数由函数值返回
	  * @param xStart - 求根区间的左端点
	  * @param xEnd - 求根区间的右端点
	  * @param dblStep - 搜索求根时采用的步长
	  * @param eps - 精度控制参数
	  * @return int 型，求得的实根的数目
	  */
	 public int getRootBisect(int nNumRoots, double[] x, double xStart, double xEnd, double dblStep, double eps)
	 {
	 	 int n,js;
	     double z,y,z1,y1,z0,y0;
	
	 	 // 根的个数清0
	     n = 0; 
	
	 	 // 从左端点开始搜索
	 	 z = xStart; 
	 	 y = func(z);
	
	 	 // 循环求解
	     while ((z<=xEnd+dblStep/2.0)&&(n!=nNumRoots))
	     { 
	 		if (Math.abs(y)<eps)
	         { 
	 			n=n+1; 
	 			x[n-1]=z;
	             z=z+dblStep/2.0; 
	 			y=func(z);
	         }
	         else
	         { 
	 			z1=z+dblStep; 
	 			y1=func(z1);
	             
	 			if (Math.abs(y1)<eps)
	             { 
	 				n=n+1; 
	 				x[n-1]=z1;
	                 z=z1+dblStep/2.0; 
	 				y=func(z);
	             }
	             else if (y*y1>0.0)
	             { 
	 				y=y1; 
	 				z=z1;
	 			}
	             else
	             { 
	 				js=0;
	                 while (js==0)
	                 { 
	 					if (Math.abs(z1-z)<eps)
	                     { 
	 						n=n+1; 
	 						x[n-1]=(z1+z)/2.0;
	                         z=z1+dblStep/2.0; y=func(z);
	                         js=1;
	                     }
	                     else
	                     { 
	 						z0=(z1+z)/2.0; 
	 						y0=func(z0);
	                         if (Math.abs(y0)<eps)
	                         { 
	 							x[n]=z0; 
	 							n=n+1; 
	 							js=1;
	                             z=z0+dblStep/2.0; 
	 							y=func(z);
	                         }
	                         else if ((y*y0)<0.0)
	                         { 
	 							z1=z0; 
	 							y1=y0;
	 						}
	                         else 
	 						{ 
	 							z=z0; 
	 							y=y0;
	 						}
	                     }
	                 }
	             }
	         }
	     }
	     
	 	// 返回实根的数目
	 	return(n);
	 }
	
	 /**
	  * 求非线性方程一个实根的牛顿法
	  * 
	  * 调用时，须覆盖计算方程左端函数f(x)及其一阶导数f'(x)值的虚函数:
	  * 	void func(double x, double[] y)
	  * 	y(0) 返回f(x)的值
	  * 	y(1) 返回f'(x)的值
	  * 
	  * @param x - 传入迭代初值（猜测解），返回在区间求得的一个实根
	  * @param nMaxIt - 递归次数
	  * @param eps - 精度控制参数
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootNewton(Real x, int nMaxIt, double eps)
	 { 
	     int l;
	     double d,p,x0,x1=0.0;
	     double[] y = new double[2];
	
         // 条件值
	     l=nMaxIt; 
	     x0=x.doubleValue();
	     func(x0,y);
	     
	 	 // 求解，控制精度
	     d=eps+1.0;
	     while ((d>=eps)&&(l!=0))
	     { 
	 		 if (y[1] == 0.0)
	 			return false;
	
	         x1=x0-y[0]/y[1];
	         func(x1,y);
	         
	 		 d=Math.abs(x1-x0); 
	 		 p=Math.abs(y[0]);
	         if (p>d) 
	 			d=p;
	         x0=x1; 
	 		 l=l-1;
	     }
	     
	 	x.setValue(x1);
	
	 	return true;
	 }
	
	 /**
	  * 求非线性方程一个实根的埃特金迭代法
	  * 
	  * 调用时，须覆盖计算方程左端函数f(x)值的虚函数
	  *       double func(double x)
	  * 
	  * @param x - 传入迭代初值（猜测解），返回在区间求得的一个实根
	  * @param nMaxIt - 递归次数
	  * @param eps - 精度控制参数
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootAitken(Real x, int nMaxIt, double eps)
	 { 
	 	int flag,l;
	    double u,v,x0;
	
	 	// 求解条件
	    l=0; 
	 	x0=x.doubleValue(); 
	 	flag=0;
	
	 	// 迭代求解
	     while ((flag==0)&&(l!=nMaxIt))
	     { 
	 		l=l+1; 
	        u=func(x0); 
	 		v=func(u);
	        if (Math.abs(u-v)<eps) 
	 		{ 
	 			x0=v; 
	 			flag=1; 
	 		}
	        else 
	 			x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
	     }
	     
	 	x.setValue(x0); 
	     
	 	// 是否在指定的迭代次数内达到求解精度
	 	return (nMaxIt > l);
	 }
	
	 /**
	  * 求非线性方程一个实根的连分式解法
	  * 
	  * 调用时，须覆盖计算方程左端函数f(x)值的虚函数
	  *       double func(double x)
	  * 
	  * @param x - 传入迭代初值（猜测解），返回在区间求得的一个实根
	  * @param eps - 精度控制参数
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootPq(Real x, double eps)
	 { 
	    int i,j,m,it=0,l;
	    double z,h,x0,q;
	    double[] a = new double[10];
	    double[] y = new double[10];
	
	 	// 求解条件
	    l=10; 
	 	q=1.0e+35; 
	 	x0=x.doubleValue(); 
	 	h=0.0;
	     
	 	// 连分式求解
	 	while (l!=0)
	    { 
	 		l=l-1; 
	 		j=0; 
	 		it=l;
	        while (j<=7)
	        { 
	 			if (j<=2) 
	 				z=x0+0.1*j;
	            else 
	 				z=h;
	 			
	 			y[j]=func(z);
	            h=z;
	            if (j==0) 
	 				a[0]=z;
	             else
	             { 
	 				m=0; 
	 				i=0;
	                while ((m==0)&&(i<=j-1))
	                { 
	 					if (Math.abs(h-a[i])+1.0==1.0) 
	 						m=1;
	                    else 
	 						h=(y[j]-y[i])/(h-a[i]);
	                      
	 					i=i+1;
	                 }
	                 a[j]=h;
	                 if (m!=0) 
	 					a[j]=q;
	                 h=0.0;
	                 for (i=j-1; i>=0; i--)
	                 { 
	 					if (Math.abs(a[i+1]+h)+1.0==1.0) 
	 						h=q;
	                    else 
	 						h=-y[i]/(a[i+1]+h);
	                 }
	                  
	 				h=h+a[0];
	             }
	              
	 			if (Math.abs(y[j])>=eps) 
	 				j=j+1;
	            else 
	 			{ 
	 				j=10; 
	 				l=0;
	 			}
	 		}
	         
	 		x0=h;
	 	}
	
	 	x.setValue(h);
	     
	 	// 是否在10阶连分式内求的实根？
	 	return (10>it);
	 }
	
	 /**
	  * 求实系数代数方程全部根的QR方法
	  * 
	  * @param n - 多项式方程的次数
	  * @param dblCoef - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的
	  *                  n+1个系数
	  * @param xr - 一维数组，长度为n，返回n个根的实部
	  * @param xi - 一维数组，长度为n，返回n个根的虚部
	  * @param nMaxIt - 迭代次数
	  * @param eps - 精度控制参数
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootQr(int n, double[] dblCoef, double[] xr, double[] xi, int nMaxIt, double eps)
	 { 
	 	// 初始化矩阵
	 	Matrix mtxQ = new Matrix();
	 	mtxQ.init(n, n);
	 	double[] q = mtxQ.getData();
	
	 	// 构造赫申伯格矩阵
	     for (int j=0; j<=n-1; j++)
	 		q[j]=-dblCoef[n-j-1]/dblCoef[n];
	
	     for (int j=n; j<=n*n-1; j++)
	 		q[j]=0.0;
	
	     for (int i=0; i<=n-2; i++)
	 		q[(i+1)*n+i]=1.0;
	
	 	// 求赫申伯格矩阵的特征值和特征向量，即为方程的解
	     if (mtxQ.computeEvHBerg(xr, xi, nMaxIt, eps))
	 		return true;
	
	 	return false;
	 }
	
	 /**
	  * 求实系数代数方程全部根的牛顿下山法
	  * 
	  * @param n - 多项式方程的次数
	  * @param dblCoef - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的
	  *                  n+1个系数
	  * @param xr - 一维数组，长度为n，返回n个根的实部
	  * @param xi - 一维数组，长度为n，返回n个根的虚部
	  * @return boolean 型，求解是否成功
	  */
	 public boolean getRootNewtonDownHill(int n, double[] dblCoef, double[] xr, double[] xi)
	 { 
	 	int m=0,i=0,jt=0,k=0,is=0,it=0;
	    double t=0,x=0,y=0,x1=0,y1=0,dx=0,dy=0,p=0,q=0,w=0,dd=0,dc=0,c=0;
	    double g=0,u=0,v=0,pq=0,g1=0,u1=0,v1=0;
	     
	 	// 初始判断
	    m=n;
	    while ((m>0)&&(Math.abs(dblCoef[m])+1.0==1.0)) 
	 		m=m-1;
	
	 	// 求解失败
	    if (m<=0)
	 		return false;
	
	    for (i=0; i<=m; i++)
	 		dblCoef[i]=dblCoef[i]/dblCoef[m];
	     
	 	for (i=0; i<=m/2; i++)
	    { 
	 		w=dblCoef[i]; 
	 		dblCoef[i]=dblCoef[m-i]; 
	 		dblCoef[m-i]=w;
	 	}
	     
	 	// 迭代求解
	 	k=m; 
	 	is=0; 
	 	w=1.0;
	    jt=1;
	    while (jt==1)
	    { 
	 		pq=Math.abs(dblCoef[k]);
	 		while (pq<1.0e-12)
	        { 
	 			xr[k-1]=0.0; 
	 			xi[k-1]=0.0; 
	 			k=k-1;
	 			if (k==1)
	            { 
	 				xr[0]=-dblCoef[1]*w/dblCoef[0]; 
	 				xi[0]=0.0;
	                 
	 				return true;
	            }
	             
	 			pq=Math.abs(dblCoef[k]);
	        }
	 	
	 		q=Math.log(pq); 
	 		q=q/(1.0*k); 
	 		q=Math.exp(q);
	        p=q; 
	 		w=w*p;
	        for (i=1; i<=k; i++)
	        {