nlequations.cs

csharp版常见数值计算源码

/*
 * 求解非线性方程组的类 NLEquations
 * 
 * 周长发编制
 */

using System;

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


		/**
		 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
		 * 
		 * @param x - 变量
		 * @param y - 函数值数组
		 */
		protected virtual void Func(double x, double[] y)
		{
		}
	
		/**
		 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
		 * 
		 * @param x - 二元函数的变量
		 * @param y - 二元函数的变量
		 * @return 函数值
		 */
		protected virtual double Func(double x, double y)
		{
			return 0;
		}
	
		/**
		 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
		 * 
		 * @param x - 二元函数的变量值数组
		 * @param y - 二元函数的变量值数组
		 * @return 函数值
		 */
		protected virtual double Func(double[] x, double[] y)
		{
			return 0;
		}
	
		/**
		 * 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
		 * 
		 * @param x - 已知变量值数组
		 * @param p - 已知函数值数组
		 */
		protected virtual 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 bool 型，求解是否成功
		 */
		public bool GetRootNewton(ref double x, int nMaxIt, double eps)
		{ 
			int l;
			double d,p,x0,x1=0.0;
			double[] y = new double[2];
	
			// 条件值
			l=nMaxIt; 
			x0=x;
			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 = x1;
	
			return true;
		}
	
		/**
		 * 求非线性方程一个实根的埃特金迭代法
		 * 
		 * 调用时，须覆盖计算方程左端函数f(x)值的虚函数
		 *       double Func(double x)
		 * 
		 * @param x - 传入迭代初值（猜测解），返回在区间求得的一个实根
		 * @param nMaxIt - 递归次数
		 * @param eps - 精度控制参数
		 * @return bool 型，求解是否成功
		 */
		public bool GetRootAitken(ref double x, int nMaxIt, double eps)
		{ 
			int flag,l;
			double u,v,x0;
	
			// 求解条件
			l=0; 
			x0=x; 
			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 = x0; 
	     
			// 是否在指定的迭代次数内达到求解精度
			return (nMaxIt > l);
		}
	
		/**
		 * 求非线性方程一个实根的连分式解法
		 * 
		 * 调用时，须覆盖计算方程左端函数f(x)值的虚函数
		 *       double Func(double x)
		 * 
		 * @param x - 传入迭代初值（猜测解），返回在区间求得的一个实根
		 * @param eps - 精度控制参数
		 * @return bool 型，求解是否成功
		 */
		public bool GetRootPq(ref double 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; 
			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 = 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 bool 型，求解是否成功
		 */
		public bool 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 bool 型，求解是否成功
		 */
		public bool GetRootNewtonDownHill(int n, double[] dblCoef, double[] xr, double[] xi)
		{ 
			int m=0,i=0,jt=0,k=0,nis=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; 
			nis=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++)
				{ 
					dblCoef[i]=dblCoef[i]/q; 
					q=q*p;
				}
	         
				x=0.0001; 
				x1=x; 
				y=0.2; 
				y1=y; 
				dx=1.0;
				g=1.0e+37;
	        
				while (true)
				{
					u=dblCoef[0]; v=0.0;
					for (i=1; i<=k; i++)
					{ 
						p=u*x1; 
						q=v*y1;
						pq=(u+v)*(x1+y1);
						u=p-q+dblCoef[i]; 
						v=pq-p-q;
					}
		         
					g1=u*u+v*v;
					if (g1>=g)
					{ 
						if (nis!=0)
						{