lequations.java

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

		// 方程组的属性，将常数矩阵赋给解矩阵
		mtxResult.setValue(mtxLEConst);
		mtxResultImag.setValue(mtxConstImag);
		double[] pDataCoef = mtxLECoef.getData();
		double[] pDataConst = mtxResult.getData();
		double[] pDataCoefImag = mtxCoefImag.getData();
		double[] pDataConstImag = mtxResultImag.getData();
		int n = getNumUnknowns();
		int m = mtxLEConst.getNumColumns();

		// 临时缓冲区，存放变换的列数
	    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++)
				{ 
					u=i*n+j;
					p=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
					if (p>d) 
					{
						d=p;
						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;
	                p=pDataCoef[u]; 
					pDataCoef[u]=pDataCoef[v]; 
					pDataCoef[v]=p;
	                p=pDataCoefImag[u]; 
					pDataCoefImag[u]=pDataCoefImag[v]; 
					pDataCoefImag[v]=p;
	            }
	            
				for (j=0;j<=m-1;j++)
	            { 
					u=k*m+j; 
					v=nIs*m+j;
	                p=pDataConst[u]; 
					pDataConst[u]=pDataConst[v]; 
					pDataConst[v]=p;
	                p=pDataConstImag[u]; 
					pDataConstImag[u]=pDataConstImag[v]; 
					pDataConstImag[v]=p;
	            }
	        }
	        
			if (pnJs[k]!=k)
	        {
				for (i=0;i<=n-1;i++)
	            { 
					u=i*n+k; 
					v=i*n+pnJs[k];
					p=pDataCoef[u]; 
					pDataCoef[u]=pDataCoef[v]; 
					pDataCoef[v]=p;
					p=pDataCoefImag[u]; 
					pDataCoefImag[u]=pDataCoefImag[v]; 
					pDataCoefImag[v]=p;
	            }
			}

	        v=k*n+k;
	        for (j=k+1;j<=n-1;j++)
	        { 
				u=k*n+j;
	            p=pDataCoef[u]*pDataCoef[v]; 
				q=-pDataCoefImag[u]*pDataCoefImag[v];
	            s=(pDataCoef[v]-pDataCoefImag[v])*(pDataCoef[u]+pDataCoefImag[u]);
	            pDataCoef[u]=(p-q)/d; 
				pDataCoefImag[u]=(s-p-q)/d;
	        }
	        
			for (j=0;j<=m-1;j++)
	        { 
				u=k*m+j;
	            p=pDataConst[u]*pDataCoef[v]; 
				q=-pDataConstImag[u]*pDataCoefImag[v];
	            s=(pDataCoef[v]-pDataCoefImag[v])*(pDataConst[u]+pDataConstImag[u]);
	            pDataConst[u]=(p-q)/d; 
				pDataConstImag[u]=(s-p-q)/d;
	        }
	        
			for (i=0;i<=n-1;i++)
	        {
				if (i!=k)
	            { 
					u=i*n+k;
					for (j=k+1;j<=n-1;j++)
	                { 
						v=k*n+j; 
						l=i*n+j;
						p=pDataCoef[u]*pDataCoef[v]; 
						q=pDataCoefImag[u]*pDataCoefImag[v];
						s=(pDataCoef[u]+pDataCoefImag[u])*(pDataCoef[v]+pDataCoefImag[v]);
						pDataCoef[l]=pDataCoef[l]-p+q;
						pDataCoefImag[l]=pDataCoefImag[l]-s+p+q;
	                }
	            
					for (j=0;j<=m-1;j++)
					{ 
						l=i*m+j; 
						v=k*m+j;
						p=pDataCoef[u]*pDataConst[v]; q=pDataCoefImag[u]*pDataConstImag[v];
						s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[v]+pDataConstImag[v]);
						pDataConst[l]=pDataConst[l]-p+q; 
						pDataConstImag[l]=pDataConstImag[l]-s+p+q;
					}
	            }
			}
	    }

		// 求解调整
		for (k=n-1;k>=0;k--)
	    {
			if (pnJs[k]!=k)
			{
				for (j=0;j<=m-1;j++)
				{ 
					u=k*m+j;
					v=pnJs[k]*m+j;
					p=pDataConst[u]; 
					pDataConst[u]=pDataConst[v]; 
					pDataConst[v]=p;
					p=pDataConstImag[u]; 
					pDataConstImag[u]=pDataConstImag[v]; 
					pDataConstImag[v]=p;
				}
			}
		}

		return true;
	}

	/**
	 * 求解三对角线方程组的追赶法
	 * 
	 * @param mtxResult - Matrix对象，返回方程组解矩阵
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetTriDiagonal(Matrix mtxResult)
	{ 
		int k,j;
	    double s;
	    
		// 将常数矩阵赋给解矩阵
		mtxResult.setValue(mtxLEConst);
		double[] pDataConst = mtxResult.getData();

		int n = getNumUnknowns();
		if (mtxLECoef.getNumRows() != n)
			return false;

		// 为系数矩阵对角线数组分配内存
		double[] pDiagData = new double[3*n-2];

		// 构造系数矩阵对角线元素数组
		k = j = 0;
		if (k == 0)
		{
			pDiagData[j++] = mtxLECoef.getElement(k, k);
			pDiagData[j++] = mtxLECoef.getElement(k, k+1);
		}
	    for (k=1; k<n-1; ++k)
		{
			pDiagData[j++] = mtxLECoef.getElement(k, k-1);
			pDiagData[j++] = mtxLECoef.getElement(k, k);
			pDiagData[j++] = mtxLECoef.getElement(k, k+1);
		}
		if (k == n-1)
		{
			pDiagData[j++] = mtxLECoef.getElement(k, k-1);
			pDiagData[j++] = mtxLECoef.getElement(k, k);
		}

		// 求解
		for (k=0;k<=n-2;k++)
	    { 
			j=3*k; 
			s=pDiagData[j];

			// 求解失败
	        if (Math.abs(s)+1.0==1.0)
			{
				return false;
			}

	        pDiagData[j+1]=pDiagData[j+1]/s;
	        pDataConst[k]=pDataConst[k]/s;
	        pDiagData[j+3]=pDiagData[j+3]-pDiagData[j+2]*pDiagData[j+1];
	        pDataConst[k+1]=pDataConst[k+1]-pDiagData[j+2]*pDataConst[k];
	    }
	    
		s=pDiagData[3*n-3];
	    if (s == 0.0)
		{
			return false;
		}
	    
		// 调整
		pDataConst[n-1]=pDataConst[n-1]/s;
	    for (k=n-2;k>=0;k--)
			pDataConst[k]=pDataConst[k]-pDiagData[3*k+1]*pDataConst[k+1];
	    
		return true;
	}

	/**
	 * 一般带型方程组的求解
	 * 
	 * @param nBandWidth - 带宽
	 * @param mtxResult - Matrix对象，返回方程组解矩阵
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetBand(int nBandWidth, Matrix mtxResult)
	{ 
		int ls,k,i,j,is=0,u,v;
	    double p,t;
	    
		// 带宽必须为奇数
		if ((nBandWidth-1)%2 != 0)
			return false;

		// 将常数矩阵赋给解矩阵
		mtxResult.setValue(mtxLEConst);
		double[] pDataConst = mtxResult.getData();

		// 方程组属性
		int m = mtxLEConst.getNumColumns();
		int n = getNumUnknowns();
		if (mtxLECoef.getNumRows() != n)
			return false;

		// 带宽数组：带型矩阵的有效数据
		double[] pBandData = new double[nBandWidth*n];

		// 半带宽
	    ls = (nBandWidth-1)/2;

		// 构造带宽数组
	    for (i=0; i<n; ++i)
		{
			j = 0;
			for (k=Math.max(0, i-ls); k<Math.max(0, i-ls)+nBandWidth; ++k)
			{
				if (k < n)
					pBandData[i*nBandWidth+j++] = mtxLECoef.getElement(i, k);
				else
					pBandData[i*nBandWidth+j++] = 0;
			}
		}

		// 求解
	    for (k=0;k<=n-2;k++)
	    { 
			p=0.0;
	        for (i=k;i<=ls;i++)
	        { 
				t=Math.abs(pBandData[i*nBandWidth]);
	            if (t>p) 
				{
					p=t; 
					is=i;
				}
	        }
	        
			if (p == 0.0)
			{
				return false;
			}

	        for (j=0;j<=m-1;j++)
	        { 
				u=k*m+j; 
				v=is*m+j;
	            t=pDataConst[u]; 
				pDataConst[u]=pDataConst[v]; 
				pDataConst[v]=t;
	        }
	        
			for (j=0;j<=nBandWidth-1;j++)
	        { 
				u=k*nBandWidth+j; 
				v=is*nBandWidth+j;
	            t=pBandData[u]; 
				pBandData[u]=pBandData[v]; 
				pBandData[v]=t;
	        }
	        
			for (j=0;j<=m-1;j++)
	        { 
				u=k*m+j; 
				pDataConst[u]=pDataConst[u]/pBandData[k*nBandWidth];
			}
	        
			for (j=1;j<=nBandWidth-1;j++)
	        { 
				u=k*nBandWidth+j; 
				pBandData[u]=pBandData[u]/pBandData[k*nBandWidth];
			}
	        
			for (i=k+1;i<=ls;i++)
	        { 
				t=pBandData[i*nBandWidth];
	            for (j=0;j<=m-1;j++)
	            { 
					u=i*m+j; 
					v=k*m+j;
	                pDataConst[u]=pDataConst[u]-t*pDataConst[v];
	            }
	            
				for (j=1;j<=nBandWidth-1;j++)
	            { 
					u=i*nBandWidth+j; 
					v=k*nBandWidth+j;
	                pBandData[u-1]=pBandData[u]-t*pBandData[v];
	            }
	            
				u=i*nBandWidth+nBandWidth-1; pBandData[u]=0.0;
	        }
	        
			if (ls!=(n-1)) 
				ls=ls+1;
	    }
	    
		p=pBandData[(n-1)*nBandWidth];
	    if (p == 0.0)
		{
			return false;
		}

	    for (j=0;j<=m-1;j++)
	    { 
			u=(n-1)*m+j; 
			pDataConst[u]=pDataConst[u]/p;
		}
	    
		ls=1;
	    for (i=n-2;i>=0;i--)
	    { 
			for (k=0;k<=m-1;k++)
	        { 
				u=i*m+k;
	            for (j=1;j<=ls;j++)
	            { 
					v=i*nBandWidth+j; 
					is=(i+j)*m+k;
	                pDataConst[u]=pDataConst[u]-pBandData[v]*pDataConst[is];
	            }
	        }
	        
			if (ls!=(nBandWidth-1)) 
				ls=ls+1;
	    }
	    
		return true;
	}

	/**
	 * 求解对称方程组的分解法
	 * 
	 * @param mtxResult - CMatrix引用对象，返回方程组解矩阵
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetDjn(Matrix mtxResult)
	{ 
		int i,j,l,k,u,v,w,k1,k2,k3;
	    double p;

		// 方程组属性，将常数矩阵赋给解矩阵
		Matrix mtxCoef = new Matrix(mtxLECoef);
		mtxResult.setValue(mtxLEConst);
		int n = mtxCoef.getNumColumns();
		int m = mtxResult.getNumColumns();
		double[] pDataCoef = mtxCoef.getData();
		double[] pDataConst = mtxResult.getData();

		// 非对称系数矩阵，不能用本方法求解
	    if (pDataCoef[0] == 0.0)
			return false;

	    for (i=1; i<=n-1; i++)
	    { 
			u=i*n; 
			pDataCoef[u]=pDataCoef[u]/pDataCoef[0];
		}
	    
		for (i=1; i<=n-2; i++)
	    { 
			u=i*n+i;
	        for (j=1; j<=i; j++)
	        { 
				v=i*n+j-1; 
				l=(j-1)*n+j-1;
	            pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v]*pDataCoef[l];
	        }
	        
			p=pDataCoef[u];
	        if (p == 0.0)
				return false;

	        for (k=i+1; k<=n-1; k++)
	        { 
				u=k*n+i;
	            for (j=1; j<=i; j++)
	            { 
					v=k*n+j-1; 
					l=i*n+j-1; 
					w=(j-1)*n+j-1;
					pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[l]*pDataCoef[w];
	            }
	            
				pDataCoef[u]=pDataCoef[u]/p;
	        }
	    }
	    
		u=n*n-1;
	    for (j=1; j<=n-1; j++)
	    { 
			v=(n-1)*n+j-1; 
			w=(j-1)*n+j-1;
	        pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v]*pDataCoef[w];
	    }
	    
		p=pDataCoef[u];
	    if (p == 0.0)
			return false;

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

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

	/**
	 * 求解对称正定方程组的平方根法
	 * 
	 * @param mtxResult - CMatrix引用对象，返回方程组解矩阵
	 * @return boolean 型，方程组求解是否成功
	 */
	public boolean getRootsetCholesky(Matrix mtxResult)
	{ 
		int i,j,k,u,v;
	    
		// 方程组属性，将常数矩阵赋给解矩阵
		Matrix mtxCoef = new Matrix(mtxLECoef);
		mtxResult.setValue(mtxLEConst);
		int n = mtxCoef.getNumColumns();
		int m = mtxResult.getNumColumns();
		double[] pDataCoef = mtxCoef.getData();
		double[] pDataConst = mtxResult.getData();
	    
		// 非对称正定系数矩阵，不能用本方法求解
		if (pDataCoef[0] <= 0.0)
			return false;

		pDataCoef[0]=Math.sqrt(pDataCoef[0]);
	    for (j=1; j<=n-1; j++) 
			pDataCoef[j]=pDataCoef[j]/pDataCoef[0];
	    
		for (i=1; i<=n-1; i++)
	    { 
			u=i*n+i;
	        for (j=1; j<=i; j++)
	        { 
				v=(j-1)*n+i;
	            pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v];
	        }
	        
			if (pDataCoef[u] <= 0.0)
				return false;

	        pDataCoef[u]=Math.sqrt(pDataCoef[u]);
	        if (i!=(n-1))
	        { 
				for (j=i+1; j<=n-1; j++)
	            { 
					v=i*n+j;
	                for (k=1; k<=i; k++)
						pDataCoef[v]=pDataCoef[v]-pDataCoef[(k-1)*n+i]*pDataCoef[(k-1)*n+j];
	                pDataCoef[v]=pDataCoef[v]/pDataCoef[u];
	            }
	        }
	    }
	    
		for (j=0; j<=m-1; j++)
	    { 
			pDataConst[j]=pDataConst[j]/pDataCoef[0];