📄 nlequations.java
字号:
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 (is!=0)
{
it=1;
Real xtmp = new Real(x);
Real ytmp = new Real(y);
Real x1tmp = new Real(x1);
Real y1tmp = new Real(y1);
Real dxtmp = new Real(dx);
Real dytmp = new Real(dy);
Real ddtmp = new Real(dd);
Real dctmp = new Real(dc);
Real ctmp = new Real(c);
Real ktmp = new Real(k);
Real istmp = new Real(is);
Real ittmp = new Real(it);
g65c(xtmp, ytmp, x1tmp, y1tmp, dxtmp, dytmp, ddtmp, dctmp, ctmp, ktmp, istmp, ittmp);
x=xtmp.doubleValue();
y=ytmp.doubleValue();
x1=x1tmp.doubleValue();
y1=y1tmp.doubleValue();
dx=dxtmp.doubleValue();
dy=dytmp.doubleValue();
dd=ddtmp.doubleValue();
dc=dctmp.doubleValue();
c=ctmp.doubleValue();
k=ktmp.intValue();
is=istmp.intValue();
it=ittmp.intValue();
if (it==0)
continue;
}
else
{
Real ttmp = new Real(t);
Real xtmp = new Real(x);
Real ytmp = new Real(y);
Real x1tmp = new Real(x1);
Real y1tmp = new Real(y1);
Real dxtmp = new Real(dx);
Real dytmp = new Real(dy);
Real ptmp = new Real(p);
Real qtmp = new Real(q);
Real ktmp = new Real(k);
Real ittmp = new Real(it);
g60c(ttmp,xtmp,ytmp,x1tmp,y1tmp,dxtmp,dytmp,ptmp,qtmp,ktmp,ittmp);
t=ttmp.doubleValue();
x=xtmp.doubleValue();
y=ytmp.doubleValue();
x1=x1tmp.doubleValue();
y1=y1tmp.doubleValue();
dx=dxtmp.doubleValue();
dy=dytmp.doubleValue();
p=ptmp.doubleValue();
q=qtmp.doubleValue();
k=ktmp.intValue();
it=ittmp.intValue();
if (t>=1.0e-03)
continue;
if (g>1.0e-18)
{
it=0;
Real xtmp1 = new Real(x);
Real ytmp1 = new Real(y);
Real x1tmp1 = new Real(x1);
Real y1tmp1 = new Real(y1);
Real dxtmp1 = new Real(dx);
Real dytmp1 = new Real(dy);
Real ddtmp1 = new Real(dd);
Real dctmp1 = new Real(dc);
Real ctmp1 = new Real(c);
Real ktmp1 = new Real(k);
Real istmp1 = new Real(is);
Real ittmp1 = new Real(it);
g65c(xtmp1, ytmp1, x1tmp1, y1tmp1, dxtmp1, dytmp1, ddtmp1, dctmp1, ctmp1, ktmp1, istmp1, ittmp1);
x=xtmp1.doubleValue();
y=ytmp1.doubleValue();
x1=x1tmp1.doubleValue();
y1=y1tmp1.doubleValue();
dx=dxtmp1.doubleValue();
dy=dytmp1.doubleValue();
dd=ddtmp1.doubleValue();
dc=dctmp1.doubleValue();
c=ctmp1.doubleValue();
k=ktmp1.intValue();
is=istmp1.intValue();
it=ittmp1.intValue();
if (it==0)
continue;
}
}
Real xtmp = new Real(x);
Real ytmp = new Real(y);
Real ptmp = new Real(p);
Real wtmp = new Real(w);
Real ktmp = new Real(k);
g90c(xr,xi,ar,ai,xtmp,ytmp,ptmp,wtmp,ktmp);
x=xtmp.doubleValue();
y=ytmp.doubleValue();
p=ptmp.doubleValue();
w=wtmp.doubleValue();
k=ktmp.intValue();
break;
}
else
{
g=g1;
x=x1;
y=y1;
is=0;
if (g<=1.0e-22)
{
Real xtmp = new Real(x);
Real ytmp = new Real(y);
Real ptmp = new Real(p);
Real wtmp = new Real(w);
Real ktmp = new Real(k);
g90c(xr,xi,ar,ai,xtmp,ytmp,ptmp,wtmp,ktmp);
x=xtmp.doubleValue();
y=ytmp.doubleValue();
p=ptmp.doubleValue();
w=wtmp.doubleValue();
k=ktmp.intValue();
}
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;
Real xtmp = new Real(x);
Real ytmp = new Real(y);
Real x1tmp = new Real(x1);
Real y1tmp = new Real(y1);
Real dxtmp = new Real(dx);
Real dytmp = new Real(dy);
Real ddtmp = new Real(dd);
Real dctmp = new Real(dc);
Real ctmp = new Real(c);
Real ktmp = new Real(k);
Real istmp = new Real(is);
Real ittmp = new Real(it);
g65c(xtmp, ytmp, x1tmp, y1tmp, dxtmp, dytmp, ddtmp, dctmp, ctmp, ktmp, istmp, ittmp);
x=xtmp.doubleValue();
y=ytmp.doubleValue();
x1=x1tmp.doubleValue();
y1=y1tmp.doubleValue();
dx=dxtmp.doubleValue();
dy=dytmp.doubleValue();
dd=ddtmp.doubleValue();
dc=dctmp.doubleValue();
c=ctmp.doubleValue();
k=ktmp.intValue();
is=istmp.intValue();
it=ittmp.intValue();
if (it==0)
continue;
Real xtmp1 = new Real(x);
Real ytmp1 = new Real(y);
Real ptmp1 = new Real(p);
Real wtmp1 = new Real(w);
Real ktmp1 = new Real(k);
g90c(xr,xi,ar,ai,xtmp1,ytmp1,ptmp1,wtmp1,ktmp1);
x=xtmp1.doubleValue();
y=ytmp1.doubleValue();
p=ptmp1.doubleValue();
w=wtmp1.doubleValue();
k=ktmp1.intValue();
}
else
{
dx=(u*u1+v*v1)/p;
dy=(u1*v-v1*u)/p;
t=1.0+4.0/k;
Real ttmp = new Real(t);
Real xtmp = new Real(x);
Real ytmp = new Real(y);
Real x1tmp = new Real(x1);
Real y1tmp = new Real(y1);
Real dxtmp = new Real(dx);
Real dytmp = new Real(dy);
Real ptmp = new Real(p);
Real qtmp = new Real(q);
Real ktmp = new Real(k);
Real ittmp = new Real(it);
g60c(ttmp,xtmp,ytmp,x1tmp,y1tmp,dxtmp,dytmp,ptmp,qtmp,ktmp,ittmp);
t=ttmp.doubleValue();
x=xtmp.doubleValue();
y=ytmp.doubleValue();
x1=x1tmp.doubleValue();
y1=y1tmp.doubleValue();
dx=dxtmp.doubleValue();
dy=dytmp.doubleValue();
p=ptmp.doubleValue();
q=qtmp.doubleValue();
k=ktmp.intValue();
it=ittmp.intValue();
if (t>=1.0e-03)
continue;
if (g>1.0e-18)
{
it=0;
Real xtmp1 = new Real(x);
Real ytmp1 = new Real(y);
Real x1tmp1 = new Real(x1);
Real y1tmp1 = new Real(y1);
Real dxtmp1 = new Real(dx);
Real dytmp1 = new Real(dy);
Real ddtmp1 = new Real(dd);
Real dctmp1 = new Real(dc);
Real ctmp1 = new Real(c);
Real ktmp1 = new Real(k);
Real istmp1 = new Real(is);
Real ittmp1 = new Real(it);
g65c(xtmp1, ytmp1, x1tmp1, y1tmp1, dxtmp1, dytmp1, ddtmp1, dctmp1, ctmp1, ktmp1, istmp1, ittmp1);
x=xtmp1.doubleValue();
y=ytmp1.doubleValue();
x1=x1tmp1.doubleValue();
y1=y1tmp1.doubleValue();
dx=dxtmp1.doubleValue();
dy=dytmp1.doubleValue();
dd=ddtmp1.doubleValue();
dc=dctmp1.doubleValue();
c=ctmp1.doubleValue();
k=ktmp1.intValue();
is=istmp1.intValue();
it=ittmp1.intValue();
if (it==0)
continue;
}
Real xtmp1 = new Real(x);
Real ytmp1 = new Real(y);
Real ptmp1 = new Real(p);
Real wtmp1 = new Real(w);
Real ktmp1 = new Real(k);
g90c(xr,xi,ar,ai,xtmp1,ytmp1,ptmp1,wtmp1,ktmp1);
x=xtmp1.doubleValue();
y=ytmp1.doubleValue();
p=ptmp1.doubleValue();
w=wtmp1.doubleValue();
k=ktmp1.intValue();
}
}
break;
}
}
if (k==1)
jt=0;
else
jt=1;
}
return true;
}
/**
* 内部函数
*/
private void g60c(Real t,Real x,Real y,Real x1,Real y1,Real dx,Real dy,Real p,
Real q,Real k,Real it)
{
it.setValue(1);
while (it.intValue()==1)
{
t.setValue(t.doubleValue()/1.67);
it.setValue(0);
x1.setValue(x.doubleValue()-t.doubleValue()*dx.doubleValue());
y1.setValue(y.doubleValue()-t.doubleValue()*dy.doubleValue());
if (k.intValue()>=30)
{
p.setValue(Math.sqrt(x1.doubleValue()*x1.doubleValue()+y1.doubleValue()*y1.doubleValue()));
q.setValue(Math.exp(75.0/k.doubleValue()));
if (p.doubleValue()>=q.doubleValue())
it.setValue(1);
}
}
}
/**
* 内部函数
*/
private void g90c(double[] xr,double[] xi,double[] ar,double[] ai,Real x,Real y,Real p,Real w,Real k)
{
int i;
for (i=1; i<=k.intValue(); i++)
{
ar[i]=ar[i]+ar[i-1]*x.doubleValue()-ai[i-1]*y.doubleValue();
ai[i]=ai[i]+ar[i-1]*y.doubleValue()+ai[i-1]*x.doubleValue();
}
xr[k.intValue()-1]=x.doubleValue()*w.doubleValue();
xi[k.intValue()-1]=y.doubleValue()*w.doubleValue();
k.setValue(k.intValue()-1);
if (k.intValue()==1)
{
p.setValue(ar[0]*ar[0]+ai[0]*ai[0]);
xr[0]=-w.doubleValue()*(ar[0]*ar[1]+ai[0]*ai[1])/p.doubleValue();
xi[0]=w.doubleValue()*(ar[1]*ai[0]-ar[0]*ai[1])/p.doubleValue();
}
}
/**
* 内部函数
*/
private void g65c(Real x,Real y,Real x1,Real y1,Real dx,Real dy,Real dd,Real dc,Real c,Real k,Real is,Real it)
{
if (it.intValue()==0)
{
is.setValue(1);
dd.setValue(Math.sqrt(dx.doubleValue()*dx.doubleValue()+dy.doubleValue()*dy.doubleValue()));
if (dd.doubleValue()>1.0)
dd.setValue(1.0);
dc.setValue(6.28/(4.5*k.doubleValue()));
c.setValue(0.0);
}
while(true)
{
c.setValue(c.doubleValue()+dc.doubleValue());
dx.setValue(dd.doubleValue()*Math.cos(c.doubleValue()));
dy.setValue(dd.doubleValue()*Math.sin(c.doubleValue()));
x1.setValue(x.doubleValue()+dx.doubleValue());
y1.setValue(y.doubleValue()+dy.doubleValue());
if (c.doubleValue()<=6.29)
{
it.setValue(0);
return;
}
dd.setValue(dd.doubleValue()/1.67);
if (dd.intValue()<=1.0e-07)
{
it.setValue(1);
return;
}
c.setValue(0.0);
}
}
/**
* 求非线性方程一个实根的蒙特卡洛法
*
* 调用时，须覆盖计算方程左端函数f(x)值的虚函数double func(double x)
*
* @param x - 传入初值（猜测解），返回求得的实根
* @param xStart - 均匀分布的端点初值
* @param nControlB - 控制参数
* @param eps - 控制精度
*/
public void getRootMonteCarlo(Real x, double xStart, int nControlB, double eps)
{
int k;
double xx,a,y,x1,y1,r;
// 求解条件
a = xStart;
k = 1;
r = 1.0;
// 初值
xx = x.doubleValue();
y = func(xx);
// 精度控制求解
while (a>=eps)
{
Real rtmp = new Real(r);
x1=rnd(rtmp);
r=rtmp.doubleValue();
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.setValue(xx);
return;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -