nlequations.cs

csharp版常见数值计算源码

							it=1;
							g65(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);
						 
							if (it==0) 
								continue;
						}
						else
						{ 
							g60(ref t,ref x,ref y,ref x1,ref y1,ref dx,ref dy,ref p,ref q,ref k,ref it);
		 				 
							if (t>=1.0e-03) 
								continue;
		                 
							if (g>1.0e-18)
							{ 
								it=0;
								g65(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);

								if (it==0) 
									continue;
							}
						}
		             
						g90(xr,xi,dblCoef,ref x,ref y,ref p,ref q,ref w,ref k);
						break;
					}
					else
					{ 
						g=g1; 
						x=x1; 
						y=y1; 
						nis=0;
						if (g<=1.0e-22)
						{
							g90(xr,xi,dblCoef,ref x,ref y,ref p,ref q,ref w,ref k);
						}
						else
						{ 
							u1=k*dblCoef[0]; 
							v1=0.0;
							for (i=2; i<=k; i++)
							{ 
								p=u1*x; 
								q=v1*y; 
								pq=(u1+v1)*(x+y);
								u1=p-q+(k-i+1)*dblCoef[i-1];
								v1=pq-p-q;
							}
		                 
							p=u1*u1+v1*v1;
							if (p<=1.0e-20)
							{ 
								it=0;
								g65(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);

								if (it==0) 
									continue;
		
								g90(xr,xi,dblCoef,ref x,ref y,ref p,ref q,ref w,ref k);
							}
							else
							{ 
								dx=(u*u1+v*v1)/p;
								dy=(u1*v-v1*u)/p;
								t=1.0+4.0/k;
								
								g60(ref t,ref x,ref y,ref x1,ref y1,ref dx,ref dy,ref p,ref q,ref k,ref it);

								if (t>=1.0e-03) 
									continue;
		
								if (g>1.0e-18)
								{ 
									it=0;

									g65(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);

									if (it==0) 
										continue;
								}
		                     
								g90(xr,xi,dblCoef,ref x,ref y,ref p,ref q,ref w,ref k);
							}
						}
						break;
					}
				}
	        
				if (k==1) 
					jt=0;
				else 
					jt=1;
			}
	     
			return true;
		}
	
	
		/**
		 * 内部函数
		 */
		private void g60(ref double t,ref double x,ref double y,ref double x1,ref double y1,ref double dx,ref double dy,ref double p,ref double q,ref int k,ref int it)
		{ 
			it=1;
			while (it==1)
			{ 
				t=t/1.67; 
				it=0;
				x1=x-(t)*(dx);
				y1=y-(t)*(dy);
				if (k>=50)
				{ 
					p=Math.Sqrt((x1)*(x1)+(y1)*(y1));
					q=Math.Exp(85.0/(k));
					if (p>=q) 
						it=1;
				}
			}
		}
	
		/**
		 * 内部函数
		 */
		private void g90(double[] xr,double[] xi,double[] dblCoef,ref double x,ref double y,ref double p,ref double q,ref double w,ref int k)
		{ 
			int i;
    
			if (Math.Abs(y)<=1.0e-06)
			{ 
				p=-(x); 
				y=0.0; 
				q=0.0;
			}
			else
			{ 
				p=-2.0*(x); 
				q=(x)*(x)+(y)*(y);
				xr[k-1]=(x)*(w);
				xi[k-1]=-(y)*(w);
				k=k-1;
			}
    
			for (i=1; i<=k; i++)
			{ 
				dblCoef[i]=dblCoef[i]-dblCoef[i-1]*(p);
				dblCoef[i+1]=dblCoef[i+1]-dblCoef[i-1]*(q);
			}
    
			xr[k-1]=(x)*(w); 
			xi[k-1]=(y)*(w);
			k=k-1;
			if (k==1)
			{ 
				xr[0]=-dblCoef[1]*(w)/dblCoef[0]; 
				xi[0]=0.0;
			}
		}
	
		/**
		 * 内部函数
		 */
		private void g65(ref double x,ref double y,ref double x1,ref double y1,ref double dx,ref double dy,ref double dd,ref double dc,ref double c,ref int k,ref int nis,ref int it)
		{ 
			if (it==0)
			{ 
				nis=1;
				dd=Math.Sqrt((dx)*(dx)+(dy)*(dy));
				if (dd>1.0) 
					dd=1.0;
				dc=6.28/(4.5*(k)); 
				c=0.0;
			}
    
			while(true)
			{ 
				c=c+(dc);
				dx=(dd)*Math.Cos(c); 
				dy=(dd)*Math.Sin(c);
				x1=x+dx; 
				y1=y+dy;
				if (c<=6.29)
				{ 
					it=0; 
					return;
				}
        
				dd=dd/1.67;
				if (dd<=1.0e-07)
				{ 
					it=1; 
					return;
				}
        
				c=0.0;
			}
		}
	
		/**
		 * 求复系数代数方程全部根的牛顿下山法
		 * 
		 * @param n - 多项式方程的次数
		 * @param ar - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的
		 *             n+1个系数的实部
		 * @param ai - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的
		 *             n+1个系数的虚部
		 * @param xr - 一维数组，长度为n，返回n个根的实部
		 * @param xi - 一维数组，长度为n，返回n个根的虚部
		 * @return bool 型，求解是否成功
		 */
		public bool GetRootNewtonDownHill(int n, double[] ar, double[] ai, 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;
			p=Math.Sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
			while ((m>0)&&(p+1.0==1.0))
			{  
				m=m-1;
				p=Math.Sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
			}
	     
			// 求解失败
			if (m<=0)
				return false;
	
			for (i=0; i<=m; i++)
			{ 
				ar[i]=ar[i]/p; 
				ai[i]=ai[i]/p;
			}
	     
			for (i=0; i<=m/2; i++)
			{ 
				w=ar[i]; 
				ar[i]=ar[m-i]; 
				ar[m-i]=w;
				w=ai[i]; 
				ai[i]=ai[m-i]; 
				ai[m-i]=w;
			}
	     
			// 迭代求解
			k=m; 
			nis=0; 
			w=1.0;
			jt=1;
			while (jt==1)
			{ 
				pq=Math.Sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
				while (pq<1.0e-12)
				{ 
					xr[k-1]=0.0; 
					xi[k-1]=0.0; 
					k=k-1;
					if (k==1)
					{ 
						p=ar[0]*ar[0]+ai[0]*ai[0];
						xr[0]=-w*(ar[0]*ar[1]+ai[0]*ai[1])/p;
						xi[0]=w*(ar[1]*ai[0]-ar[0]*ai[1])/p;
	                 
						return true;
					}
	             
					pq=Math.Sqrt(ar[k]*ar[k]+ai[k]*ai[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++)
				{ 
					ar[i]=ar[i]/q; 
					ai[i]=ai[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=ar[0]; 
					v=ai[0];
					for (i=1; i<=k; i++)
					{ 
						p=u*x1; 
						q=v*y1;
						pq=(u+v)*(x1+y1);
						u=p-q+ar[i]; 
						v=pq-p-q+ai[i];
					}
		         
					g1=u*u+v*v;
					if (g1>=g)
					{ 
						if (nis!=0)
						{ 
							it=1;
							g65c(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);
							if (it==0) 
								continue;
						}
						else
						{ 
							g60c(ref t,ref x,ref y,ref x1,ref y1,ref dx,ref dy,ref p,ref q,ref k,ref it);
							if (t>=1.0e-03) 
								continue;
		                 
							if (g>1.0e-18)
							{ 
								it=0;
								g65c(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);
								if (it==0) 
									continue;
							}
						}
		             
						g90c(xr,xi,ar,ai,ref x,ref y,ref p,ref w,ref k);
						break;
					}
					else
					{ 
						g=g1; 
						x=x1; 
						y=y1; 
						nis=0;
						if (g<=1.0e-22)
						{
							g90c(xr,xi,ar,ai,ref x,ref y,ref p,ref w,ref k);
						}
						else
						{ 
							u1=k*ar[0]; 
							v1=ai[0];
							for (i=2; i<=k; i++)
							{ 
								p=u1*x; 
								q=v1*y; 
								pq=(u1+v1)*(x+y);
								u1=p-q+(k-i+1)*ar[i-1];
								v1=pq-p-q+(k-i+1)*ai[i-1];
							}
		                 
							p=u1*u1+v1*v1;
							if (p<=1.0e-20)
							{ 
								it=0;
								g65c(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);
								if (it==0) 
									continue;
		                     
								g90c(xr,xi,ar,ai,ref x,ref y,ref p,ref w,ref k);
							}
							else
							{ 
								dx=(u*u1+v*v1)/p;
								dy=(u1*v-v1*u)/p;
								t=1.0+4.0/k;
								g60c(ref t,ref x,ref y,ref x1,ref y1,ref dx,ref dy,ref p,ref q,ref k,ref it);
								if (t>=1.0e-03) 
									continue;
		                     
								if (g>1.0e-18)
								{ 
									it=0;
									g65c(ref x, ref y, ref x1, ref y1, ref dx, ref dy, ref dd, ref dc, ref c, ref k, ref nis, ref it);
									if (it==0) 
										continue;
								}
		                     
								g90c(xr,xi,ar,ai,ref x,ref y,ref p,ref w,ref k);
							}
						}
						break;
					}
				}
	         
				if (k==1) 
					jt=0;
				else 
					jt=1;
			}
	     
			return true;
		}
	
		/**
		 * 内部函数
		 */
		private void g60c(ref double t,ref double x,ref double y,ref double x1,ref double y1,ref double dx,ref double dy,ref double p,ref double q,ref int k,ref int it)
		{ 
			it=1;
			while (it==1)
			{ 
				t=t/1.67; 
				it=0;
				x1=x-t*(dx);
				y1=y-t*(dy);
				if (k>=30)
				{ 
					p=Math.Sqrt(x1*(x1)+y1*(y1));
					q=Math.Exp(75.0/(k));
					if (p>=q) 
						it=1;
				}
			}
		}

	
		/**
		 * 内部函数
		 */
		private void g90c(double[] xr,double[] xi,double[] ar,double[] ai,ref double x,ref double y,ref double p,ref double w,ref int k)
		{ 
			int i;
			for (i=1; i<=k; i++)
			{ 
				ar[i]=ar[i]+ar[i-1]*(x)-ai[i-1]*(y);
				ai[i]=ai[i]+ar[i-1]*(y)+ai[i-1]*(x);
			}
    
			xr[k-1]=x*(w); 
			xi[k-1]=y*(w);
			k=k-1;
			if (k==1)
			{ 
				p=ar[0]*ar[0]+ai[0]*ai[0];
				xr[0]=-w*(ar[0]*ar[1]+ai[0]*ai[1])/(p);
				xi[0]=w*(ar[1]*ai[0]-ar[0]*ai[1])/(p);
			}
		}
	
		/**
		 * 内部函数
		 */
		private void g65c(ref double x,ref double y,ref double x1,ref double y1,ref double dx,ref double dy,ref double dd,ref double dc,ref double c,ref int k,ref int nis,ref int it)
		{ 
			if (it==0)
			{ 
				nis=1;
				dd=Math.Sqrt(dx*(dx)+dy*(dy));
				if (dd>1.0) 
					dd=1.0;
				dc=6.28/(4.5*(k)); 
				c=0.0;
			}
    
			while(true)
			{ 
				c=c+dc;
				dx=dd*Math.Cos(c); 
				dy=dd*Math.Sin(c);
				x1=x+dx; 
				y1=y+dy;
				if (c<=6.29)
				{ 
					it=0; 
					return;
				}
        
				dd=dd/1.67;
				if (dd<=1.0e-07)
				{ 
					it=1; 
					return;
				}
        
				c=0.0;
			}
		}
	
		/**
		 * 求非线性方程一个实根的蒙特卡洛法
		 * 
		 * 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
		 * 
		 * @param x - 传入初值（猜测解），返回求得的实根
		 * @param xStart - 均匀分布的端点初值
		 * @param nControlB - 控制参数
		 * @param eps - 控制精度
		 */
		public void GetRootMonteCarlo(ref double x, double xStart, int nControlB, double eps)
		{ 
			int k;
			double xx,a,y,x1,y1,r;