lequations.java

实现用于工程计算的的数学方法

	        for (i=1; i<=n-1; i++)
	        { 
				u=i*n+i; 
				v=i*m+j;
	            for (k=1; k<=i; k++)
					pDataConst[v]=pDataConst[v]-pDataCoef[(k-1)*n+i]*pDataConst[(k-1)*m+j];
	            pDataConst[v]=pDataConst[v]/pDataCoef[u];
	        }
	    }
	    
		for (j=0; j<=m-1; j++)
	    { 
			u=(n-1)*m+j;
	        pDataConst[u]=pDataConst[u]/pDataCoef[n*n-1];
	        for (k=n-1; k>=1; k--)
	        { 
				u=(k-1)*m+j;
	            for (i=k; i<=n-1; i++)
	            { 
					v=(k-1)*n+i;
	                pDataConst[u]=pDataConst[u]-pDataCoef[v]*pDataConst[i*m+j];
	            }
	            
				v=(k-1)*n+k-1;
	            pDataConst[u]=pDataConst[u]/pDataCoef[v];
	        }
	    }
	    
		return true;
	}

	/**
	 * 求解大型稀疏方程组的全选主元高斯－约去消去法
	 * 
	 * @param mtxResult - CMatrix引用对象，返回方程组解矩阵
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetGgje(Matrix mtxResult)
	{ 
		int i,j,k,nIs=0,u,v;
	    double d,t;
	    
		// 方程组属性，将常数矩阵赋给解矩阵
		Matrix mtxCoef = new Matrix(mtxLECoef);
		mtxResult.setValue(mtxLEConst);
		int n = mtxCoef.getNumColumns();
		double[] pDataCoef = mtxCoef.getData();
		double[] pDataConst = mtxResult.getData();

		// 临时缓冲区，存放变换的列数
		int[] pnJs = new int[n];

		// 消元
	    for (k=0; k<=n-1; k++)
	    { 
			d=0.0;
	        for (i=k; i<=n-1; i++)
			{
				for (j=k; j<=n-1; j++)
				{ 
					t=Math.abs(pDataCoef[i*n+j]);
					if (t>d) 
					{
						d=t; 
						pnJs[k]=j; 
						nIs=i;
					}
				}
			}

	        if (d == 0.0)
			{
				return false;
			}

	        if (nIs!=k)
	        { 
				for (j=k; j<=n-1; j++)
	            { 
					u=k*n+j; 
					v=nIs*n+j;
	                t=pDataCoef[u]; 
					pDataCoef[u]=pDataCoef[v]; 
					pDataCoef[v]=t;
	            }
	            
				t=pDataConst[k]; 
				pDataConst[k]=pDataConst[nIs]; 
				pDataConst[nIs]=t;
	        }
	        
			if (pnJs[k]!=k)
	        {
				for (i=0; i<=n-1; i++)
	            { 
					u=i*n+k; 
					v=i*n+pnJs[k];
					t=pDataCoef[u]; 
					pDataCoef[u]=pDataCoef[v]; 
					pDataCoef[v]=t;
	            }
			}

	        t=pDataCoef[k*n+k];
	        for (j=k+1; j<=n-1; j++)
	        { 
				u=k*n+j;
	            if (pDataCoef[u]!=0.0) 
					pDataCoef[u]=pDataCoef[u]/t;
	        }
	        
			pDataConst[k]=pDataConst[k]/t;
	        for (j=k+1; j<=n-1; j++)
	        { 
				u=k*n+j;
	            if (pDataCoef[u]!=0.0)
	            { 
					for (i=0; i<=n-1; i++)
	                { 
						v=i*n+k;
	                    if ((i!=k)&&(pDataCoef[v]!=0.0))
	                    { 
							nIs=i*n+j;
	                        pDataCoef[nIs]=pDataCoef[nIs]-pDataCoef[v]*pDataCoef[u];
	                    }
	                }
	            }
	        }
	        
			for (i=0; i<=n-1; i++)
	        { 
				u=i*n+k;
	            if ((i!=k)&&(pDataCoef[u]!=0.0))
					pDataConst[i]=pDataConst[i]-pDataCoef[u]*pDataConst[k];
	        }
	    }
	    
		// 调整
		for (k=n-1; k>=0; k--)
	    {
			if (k!=pnJs[k])
	        { 
				t=pDataConst[k]; 
				pDataConst[k]=pDataConst[pnJs[k]]; 
				pDataConst[pnJs[k]]=t;
			}
		}

		return true;
	}

	/**
	 * 求解托伯利兹方程组的列文逊方法
	 * 
	 * @param mtxResult - CMatrix引用对象，返回方程组解矩阵
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetTlvs(Matrix mtxResult)
	{ 
		int i,j,k;
	    double a,beta,q,c,h;

		// 未知数个数
	    int n = mtxLECoef.getNumColumns();

		// 初始化解解向量
		mtxResult.init(n, 1);
		double[] x = mtxResult.getData();

		// 常数数组
		double[] pDataConst = mtxLEConst.getData();

		// 建立T数组
		double[] t = new double[n];

		// 构造T数组
	    for (i=0; i<n; ++i)
			t[i] = mtxLECoef.getElement(0, i);

		// 临时数组
		double[] s = new double[n];
		double[] y = new double[n];

		// 非托伯利兹方程组，不能用本方法求解
	    a=t[0];
	    if (a == 0.0)
	    { 
			return false;
		}

		// 列文逊方法求解
	    y[0]=1.0; 
		x[0]=pDataConst[0]/a;
	    for (k=1; k<=n-1; k++)
	    { 
			beta=0.0; 
			q=0.0;
	        for (j=0; j<=k-1; j++)
	        { 
				beta=beta+y[j]*t[j+1];
	            q=q+x[j]*t[k-j];
	        }
	        
			if (a == 0.0)
			{ 
				return false;
			}

	        c=-beta/a; 
			s[0]=c*y[k-1]; 
			y[k]=y[k-1];
	        if (k!=1)
	        {
				for (i=1; i<=k-1; i++)
					s[i]=y[i-1]+c*y[k-i-1];
			}

	        a=a+c*beta;
			if (a == 0.0)
			{ 
				return false;
			}

	        h=(pDataConst[k]-q)/a;
	        for (i=0; i<=k-1; i++)
	        { 
				x[i]=x[i]+h*s[i]; 
				y[i]=s[i];
			}
	        
			x[k]=h*y[k];
	    }
	    
	    return true;
	}

	/**
	 * 高斯－赛德尔迭代法
	 *  
	 * @param mtxResult - CMatrix引用对象，返回方程组解矩阵
	 * @param eps - 控制精度
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetGaussSeidel(Matrix mtxResult, double eps)
	{ 
		int i,j,u,v;
	    double p,t,s,q;

		// 未知数个数
	    int n = mtxLECoef.getNumColumns();

		// 初始化解向量
		mtxResult.init(n, 1);
		double[] x = mtxResult.getData();

		// 系数与常数
		double[] pDataCoef = mtxLECoef.getData();
		double[] pDataConst = mtxLEConst.getData();
	    
		// 求解
		for (i=0; i<=n-1; i++)
	    { 
			u=i*n+i; 
			p=0.0; 
			x[i]=0.0;
	        for (j=0; j<=n-1; j++)
	        {
				if (i!=j)
	            { 
					v=i*n+j; 
					p=p+Math.abs(pDataCoef[v]);
				}
			}

	        if (p>=Math.abs(pDataCoef[u]))
				return false;
	    }

		// 精度控制
	    p=eps+1.0;
	    while (p>=eps)
	    { 
			p=0.0;
	        for (i=0; i<=n-1; i++)
	        { 
				t=x[i]; 
				s=0.0;
	            for (j=0; j<=n-1; j++)
					if (j!=i) 
						s=s+pDataCoef[i*n+j]*x[j];

	            x[i]=(pDataConst[i]-s)/pDataCoef[i*n+i];
	            q=Math.abs(x[i]-t)/(1.0+Math.abs(x[i]));
	            if (q>p) 
					p=q;
	        }
	    }
	    
		return true;
	}

	/**
	 * 求解对称正定方程组的共轭梯度法
	 * 
	 * @param mtxResult - CMatrix引用对象，返回方程组解矩阵
	 * @param eps - 控制精度
	 */
	public void getRootsetGrad(Matrix mtxResult, double eps)
	{ 
		int i,k;
	    double alpha,beta,d,e;

		// 未知数个数
	    int n = getNumUnknowns();

		// 初始化解向量
		mtxResult.init(n, 1);
		double[] x = mtxResult.getData();

		// 构造临时矩阵
		Matrix mtxP = new Matrix(n, 1);
		double[] p = mtxP.getData();

		double[] pDataCoef = mtxLECoef.getData();
		double[] pDataConst = mtxLEConst.getData();

	    double[] r = new double[n];

	    for (i=0; i<=n-1; i++)
	    { 
			x[i]=0.0; 
			p[i]=pDataConst[i]; 
			r[i]=pDataConst[i]; 
		}
	    
		i=0;
	    while (i<=n-1)
	    { 
			Matrix mtxS = mtxLECoef.multiply(mtxP);
			double[] s = mtxS.getData();
	        
			d=0.0; 
			e=0.0;
	        for (k=0; k<=n-1; k++)
	        { 
				d=d+p[k]*pDataConst[k]; 
				e=e+p[k]*s[k]; 
			}
	        
			alpha=d/e;
	        for (k=0; k<=n-1; k++)
				x[k]=x[k]+alpha*p[k];
	        
			Matrix mtxQ = mtxLECoef.multiply(mtxResult);
		    double[] q = mtxQ.getData();
	        
			d=0.0;
	        for (k=0; k<=n-1; k++)
	        { 
				r[k]=pDataConst[k]-q[k]; 
				d=d+r[k]*s[k]; 
			}
	        
			beta=d/e; d=0.0;
	        for (k=0; k<=n-1; k++) 
				d=d+r[k]*r[k];
	        
			// 满足精度，求解结束
			d=Math.sqrt(d);
	        if (d<eps)
				break;

	        for (k=0; k<=n-1; k++)
				p[k]=r[k]-beta*p[k];
	        
			i=i+1;
	    }
	}

	/**
	 * 求解线性最小二乘问题的豪斯荷尔德变换法
	 * 
	 * @param mtxResult - Matrix对象，返回方程组解矩阵
	 * @param mtxQ - Matrix对象，返回豪斯荷尔德变换的Q矩阵
	 * @param mtxR - Matrix对象，返回豪斯荷尔德变换的R矩阵
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetMqr(Matrix mtxResult, Matrix mtxQ, Matrix mtxR)
	{ 
		int i,j;
	    double d;

		// 方程组的方程数和未知数个数
	    int m = mtxLECoef.getNumRows();
	    int n = mtxLECoef.getNumColumns();
		// 奇异方程组
		if (m < n)
			return false;

		// 将解向量初始化为常数向量
		mtxResult.setValue(mtxLEConst);
		double[] pDataConst = mtxResult.getData();

		// 构造临时矩阵，用于QR分解
		mtxR.setValue(mtxLECoef);
		double[] pDataCoef = mtxR.getData();

		// QR分解
	    if (! mtxR.splitQR(mtxQ))
			return false;

		// 临时缓冲区
		double[] c = new double[n];
		double[] q = mtxQ.getData();

		// 求解
	    for (i=0; i<=n-1; i++)
	    { 
			d=0.0;
	        for (j=0; j<=m-1; j++)
				d=d+q[j*m+i]*pDataConst[j];
	    
			c[i]=d;
	    }
	    
		pDataConst[n-1]=c[n-1]/pDataCoef[n*n-1];
	    for (i=n-2; i>=0; i--)
	    { 
			d=0.0;
	        for (j=i+1; j<=n-1; j++)
				d=d+pDataCoef[i*n+j]*pDataConst[j];
	        
			pDataConst[i]=(c[i]-d)/pDataCoef[i*n+i];
	    }
	    
		return true;
	}

	/**
	 * 求解线性最小二乘问题的广义逆法
	 * 
	 * @param mtxResult - Matrix对象，返回方程组解矩阵
	 * @param mtxAP - Matrix对象，返回系数矩阵的广义逆矩阵
	 * @param mtxU - Matrix对象，返回U矩阵
	 * @param mtxV - Matrix对象，返回V矩阵
	 * @param eps - 控制精度
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetGinv(Matrix mtxResult, Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps)
	{ 
		int i,j;
	    
		// 方程个数和未知数个数
	    int m = mtxLECoef.getNumRows();
	    int n = mtxLECoef.getNumColumns();

		// 初始化解向量
		mtxResult.init(n, 1);

		double[] pDataConst = mtxLEConst.getData();
		double[] x = mtxResult.getData();

		// 临时矩阵
		Matrix mtxA = new Matrix(mtxLECoef);

		// 求广义逆矩阵
	    if (! mtxA.invertUV(mtxAP, mtxU, mtxV, eps))
			return false;

		double[] pAPData = mtxAP.getData();

		// 求解
	    for (i=0; i<=n-1; i++)
	    { 
			x[i]=0.0;
	        for (j=0; j<=m-1; j++)
				x[i]=x[i]+pAPData[i*m+j]*pDataConst[j];
	    }
	    
		return true;
	}

	/**
	 * 
	 * @param mtxResult - Matrix对象，返回方程组解矩阵
	 * @param nMaxIt - 叠加次数
	 * @param eps - 控制精度
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetMorbid(Matrix mtxResult, int nMaxIt /*= 60*/, double eps)
	{ 
		int i, k;
	    double q, qq;
	    
		// 方程的阶数
	    int n = getNumUnknowns();

		// 设定迭代次数, 缺省为60
		i = nMaxIt;
	    
		// 用全选主元高斯消元法求解
		LEquations leqs = new LEquations(mtxLECoef, mtxLEConst);
		if (! leqs.getRootsetGauss(mtxResult))
			return false;
		double[] x = mtxResult.getData();
	    
		q=1.0+eps;
	    while (q>=eps)
	    { 
			// 迭代次数已达最大值，仍为求得结果，求解失败
			if (i==0)
				return false;
	        
			// 迭代次数减1
			i=i-1;
			
			// 矩阵运算
			Matrix mtxE = mtxLECoef.multiply(mtxResult);
			Matrix mtxR = mtxLEConst.subtract(mtxE);

			// 用全选主元高斯消元法求解
			leqs = new LEquations(mtxLECoef, mtxR);
			Matrix mtxRR = new Matrix();
			if (! leqs.getRootsetGauss(mtxRR))
				return false;

			double[] r = mtxRR.getData();
	        
			q=0.0;
	        for ( k=0; k<=n-1; k++)
	        { 
				qq=Math.abs(r[k])/(1.0+Math.abs(x[k]+r[k]));
	            if (qq>q) 
					q=qq;
	        }
	        
			for ( k=0; k<=n-1; k++) 
				x[k]=x[k]+r[k];

	    }
	    
		// 求解成功
		return true;
	}
}