?? nlequations.cs
字號:
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;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -