lequations.cs

csharp版常见数值计算源码

			for (j=0; j<=m-1; j++)
			{ 
				pDataConst[j]=pDataConst[j]/pDataCoef[0];
				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 bool 型，方程组求解是否成功
		 */
		public bool 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 bool 型，方程组求解是否成功
		 */
		public bool 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 bool 型，方程组求解是否成功
		 */
		public bool 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 bool 型，方程组求解是否成功
		 */
		public bool 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 bool 型，方程组求解是否成功
		 */
		public bool 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 bool 型，方程组求解是否成功
		 */
		public bool 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;
		}
	}
}