nlequations.cs

csharp版常见数值计算源码

			// 求解条件
			a = xStart; 
			k = 1; 
			r = 1.0;
	
			// 初值
			xx = x; 
			y = Func(xx);
	
			// 精度控制求解
			while (a>=eps)
			{ 
				x1=rnd(ref r);
	
				x1=-a+2.0*a*x1;
				x1=xx+x1; 
				y1=Func(x1);
	         
				k=k+1;
				if (Math.Abs(y1)>=Math.Abs(y))
				{ 
					if (k>nControlB) 
					{ 
						k=1; 
						a=a/2.0; 
					}
				}
				else
				{ 
					k=1; 
					xx=x1; 
					y=y1;
					if (Math.Abs(y)<eps)
					{ 
						x = xx; 
						return; 
					}
				}
			}
	     
			x = xx; 
		}
	
		/**
		 * 内部函数
		 */
		private double rnd(ref double r)
		{
			int m;
			double s,u,v,p;
    
			s=65536.0; 
			u=2053.0; 
			v=13849.0;
			m=(int)(r/s); 
			r=r-m*s;
			r=u*r+v; 
			m=(int)(r/s);
			r=r-m*s; p=r/s;
    
			return(p);
		}
	
		/**
		 * 求实函数或复函数方程一个复根的蒙特卡洛法
		 * 
		 * 调用时，须覆盖计算方程左端函数的模值||f(x, y)||的虚函数
		 *          double Func(double x, double y)
		 * 
		 * @param x - 传入初值（猜测解）的实部，返回求得的根的实部
		 * @param y - 传入初值（猜测解）的虚部，返回求得的根的虚部
		 * @param xStart - 均匀分布的端点初值
		 * @param nControlB - 控制参数
		 * @param eps - 控制精度
		 */
		public void GetRootMonteCarlo(ref double x, ref double y, double xStart, int nControlB, double eps)
		{ 
			int k;
			double xx,yy,a,r,z,x1,y1,z1;
	
			// 求解条件与初值
			a=xStart; 
			k=1; 
			r=1.0; 
			xx=x; 
			yy=y;
			z=Func(xx,yy);
	     
			// 精度控制求解
			while (a>=eps)
			{ 
				x1=-a+2.0*a*rnd(ref r); 

				x1=xx+x1; 
	 		
				y1=-a+2.0*a*rnd(ref r);
	        
				y1=yy+y1;
	         
				z1=Func(x1,y1);
	         
				k=k+1;
				if (z1>=z)
				{ 
					if (k>nControlB) 
					{ 
						k=1; 
						a=a/2.0; 
					}
				}
				else
				{ 
					k=1; 
					xx=x1; 
					yy=y1;  
					z=z1;
					if (z<eps)
					{ 
						x = xx; 
						y = yy; 
						return; 
					}
				}
			}
	     
			x = xx; 
			y = yy; 
		}
	
		/**
		 * 求非线性方程组一组实根的梯度法
		 * 
		 * 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
		 *          double Func(double x[], double[] y)
		 * 
		 * @param n - 方程的个数，也是未知数的个数
		 * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
		 *            返回时存放方程组的一组实根
		 * @param nMaxIt - 迭代次数
		 * @param eps - 控制精度
		 * @return bool 型，求解是否成功
		 */
		public bool GetRootsetGrad(int n, double[] x, int nMaxIt, double eps)
		{ 
			int l,j;
			double f,d,s;
			double[] y = new double[n];
	
			l=nMaxIt;
			f=Func(x,y);
	
			// 控制精度，迭代求解
			while (f>=eps)
			{ 
				l=l-1;
				if (l==0) 
				{ 
					return true;
				}
	         
				d=0.0;
				for (j=0; j<=n-1; j++) 
					d=d+y[j]*y[j];
				if (d+1.0==1.0) 
				{ 
					return false;
				}
	         
				s=f/d;
				for (j=0; j<=n-1; j++) 
					x[j]=x[j]-s*y[j];
	         
				f=Func(x,y);
			}
	     
			// 是否在有效迭代次数内达到精度
			return (nMaxIt>l);
		}
	
		/**
		 * 求非线性方程组一组实根的拟牛顿法
		 * 
		 * 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
		 *          double Func(double[] x, double[] y)
		 * 
		 * @param n - 方程的个数，也是未知数的个数
		 * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
		 *            返回时存放方程组的一组实根
		 * @param t - 控制h大小的变量，0<t<1
		 * @param h - 增量初值
		 * @param nMaxIt - 迭代次数
		 * @param eps - 控制精度
		 * @return bool 型，求解是否成功
		 */
		public bool GetRootsetNewton(int n, double[] x, double t, double h, int nMaxIt, double eps)
		{ 
			int i,j,l;
			double am,z,beta,d;
	
			double[] y = new double[n];
	 	
			// 构造矩阵
			Matrix mtxCoef = new Matrix(n, n);
			Matrix mtxConst = new Matrix(n, 1);
			double[] a = mtxCoef.GetData();
			double[] b = mtxConst.GetData();
	
			// 迭代求解
			l=nMaxIt; 
			am=1.0+eps;
			while (am>=eps)
			{ 
				Func(x,b);
	         
				am=0.0;
				for (i=0; i<=n-1; i++)
				{ 
					z=Math.Abs(b[i]);
					if (z>am) 
						am=z;
				}
	         
				if (am>=eps)
				{ 
					l=l-1;
					if (l==0)
					{ 
						return false;
					}
	             
					for (j=0; j<=n-1; j++)
					{ 
						z=x[j]; 
						x[j]=x[j]+h;
	                 
						Func(x,y);
	                 
						for (i=0; i<=n-1; i++) 
							a[i*n+j]=y[i];
	                 
						x[j]=z;
					}
	             
					// 调用全选主元高斯消元法
					LEquations leqs = new LEquations(mtxCoef, mtxConst);
					Matrix mtxResult = new Matrix();
					if (! leqs.GetRootsetGauss(mtxResult))
					{
						return false;
					}
	
					mtxConst.SetValue(mtxResult);
					b = mtxConst.GetData();
	
					beta=1.0;
					for (i=0; i<=n-1; i++) 
						beta=beta-b[i];
	             
					if (beta == 0.0)
					{ 
						return false;
					}
	             
					d=h/beta;
					for (i=0; i<=n-1; i++) 
						x[i]=x[i]-d*b[i];
	             
					h=t*h;
				}
			}
	     
			// 是否在有效迭代次数内达到精度
			return (nMaxIt>l);
		}
	
		/**
		 * 求非线性方程组最小二乘解的广义逆法
		 * 
		 * 调用时，1. 须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
		 *              double Func(double[] x, double[] y)
		 *        2. 须覆盖计算雅可比矩阵函数的虚函数
		 *              double FuncMJ(double[] x, double[] y)
		 * 
		 * @param m - 方程的个数
		 * @param n - 未知数的个数
		 * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，要求不全为0，
		 * 			返回时存放方程组的最小二乘解，当m=n时，即是非线性方程组的解
		 * @param eps1 - 最小二乘解的精度控制精度
		 * @param eps2 - 奇异值分解的精度控制精度
		 * @return bool 型，求解是否成功
		 */
		public bool GetRootsetGinv(int m, int n, double[] x, double eps1, double eps2)
		{ 
			int i,j,k,l,kk,jt;
			double alpha,z=0,h2,y1,y2,y3,y0,h1;
			double[] p,d,dx;
	     
			double[] y = new double[10];
			double[] b = new double[10];
	     
			// 控制参数
			int ka = Math.Max(m, n)+1;
			double[] w = new double[ka];
	     
			// 设定迭代次数为60，迭代求解
			l=60; 
			alpha=1.0;
			while (l>0)
			{ 
				Matrix mtxP = new Matrix(m, n);
				Matrix mtxD = new Matrix(m, 1);
				p = mtxP.GetData();
				d = mtxD.GetData();
	
				Func(x,d);
				FuncMJ(x,p);
	
				// 构造线性方程组
				LEquations leqs = new LEquations(mtxP, mtxD);
				// 临时矩阵
				Matrix mtxAP = new Matrix();
				Matrix mtxU = new Matrix();
				Matrix mtxV = new Matrix();
				// 解矩阵
				Matrix mtxDX = new Matrix();
				// 基于广义逆的最小二乘解
				if (! leqs.GetRootsetGinv(mtxDX, mtxAP, mtxU, mtxV, eps2))
				{ 
					return false;
				}
	         
				dx = mtxDX.GetData();
	
				j=0; 
				jt=1; 
				h2=0.0;
				while (jt==1)
				{ 
					jt=0;
					if (j<=2) 
						z=alpha+0.01*j;
					else 
						z=h2;
	             
					for (i=0; i<=n-1; i++) 
						w[i]=x[i]-z*dx[i];
	             
					Func(w,d);
	             
					y1=0.0;
					for (i=0; i<=m-1; i++) 
						y1=y1+d[i]*d[i];
					for (i=0; i<=n-1; i++)
						w[i]=x[i]-(z+0.00001)*dx[i];
	             
					Func(w,d);
	             
					y2=0.0;
					for (i=0; i<=m-1; i++) 
						y2=y2+d[i]*d[i];
	             
					y0=(y2-y1)/0.00001;
	             
					if (Math.Abs(y0)>1.0e-10)
					{ 
						h1=y0; h2=z;
						if (j==0) 
						{ 
							y[0]=h1; 
							b[0]=h2;
						}
						else
						{ 
							y[j]=h1; 
							kk=0; 
							k=0;
							while ((kk==0)&&(k<=j-1))
							{ 
								y3=h2-b[k];
								if (Math.Abs(y3)+1.0==1.0) 
									kk=1;
								else 
									h2=(h1-y[k])/y3;
	                         
								k=k+1;
							}
	                     
							b[j]=h2;
							if (kk!=0) 
								b[j]=1.0e+35;
	                     
							h2=0.0;
							for (k=j-1; k>=0; k--)
								h2=-y[k]/(b[k+1]+h2);
	                     
							h2=h2+b[0];
						}
	                 
						j=j+1;
						if (j<=7) 
							jt=1;
						else 
							z=h2;
					}
				}
	         
				alpha=z; 
				y1=0.0; 
				y2=0.0;
				for (i=0; i<=n-1; i++)
				{ 
					dx[i]=-alpha*dx[i];
					x[i]=x[i]+dx[i];
					y1=y1+Math.Abs(dx[i]);
					y2=y2+Math.Abs(x[i]);
				}
	         
				// 求解成功
				if (y1<eps1*y2)
				{ 
					return true;
				}
	         
				l=l-1;
			}
	     
			// 求解失败
			return false;
		}
	
		/**
		 * 求非线性方程组一组实根的蒙特卡洛法
		 * 
		 * 调用时，须覆盖计算方程左端模函数值||F||的虚函数
		 *          double Func(int n, double[] x)
		 * 其返回值为Sqr(f1*f1 + f2*f2 + … + fn*fn)
		 * 
		 * @param n - 方程的个数，也是未知数的个数
		 * @param x - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
		 *            返回时存放方程组的一组实根
		 * @param xStart - 均匀分布的端点初值
		 * @param nControlB - 控制参数
		 * @param eps - 控制精度
		 */
		public void GetRootsetMonteCarlo(int n, double[] x, double xStart, int nControlB, double eps)
		{ 
			int k,i;
			double a,r,z,z1;
	
			double[] y = new double[n];
	     
			// 初值
			a=xStart; 
			k=1; 
			r=1.0; 
	
			z = Func(x);
	
			// 用精度控制迭代求解
			while (a>=eps)
			{ 
				for (i=0; i<=n-1; i++)
				{
					y[i]=-a+2.0*a*rnd(ref r)+x[i];
				}
				z1=Func(y);
	         
				k=k+1;
				if (z1>=z)
				{ 
					if (k>nControlB) 
					{ 
						k=1; 
						a=a/2.0; 
					}
				}
				else
				{ 
					k=1; 
					for (i=0; i<=n-1; i++) 
						x[i]=y[i];
	             
					// 求解成功
					z=z1;
					if (z<eps)  
					{
						return;
					}
				}
			}
		}
	}
}