【Java开源代码栏目提醒】:网学会员--在 Java开源代码编辑为广大网友搜集整理了:d11r8F.java绩等信息,祝愿广大网友取得需要的信息,参考学习。
public class d11r8F
{
void dfunc(double x[], double df[])
{
df[1] = bessj1(x[1] - 0.5) * bessj0(x[2] - 0.5) * bessj0(x[3] - 0.5);
df[2] = bessj0(x[1] - 0.5) * bessj1(x[2] - 0.5) * bessj0(x[3] - 0.5);
df[3] = bessj0(x[1] - 0.5) * bessj0(x[2] - 0.5) * bessj1(x[3] - 0.5);
}
double func2(double x[],int n)
{
return 1.0 - bessj0(x[1] - 0.5) * bessj0(x[2] - 0.5) * bessj0(x[3] - 0.5);
}
double func(double x)
{
return f1dim(x);
}
double f1dim(double x)
{
double xt[] = new double[51];
for (int j = 1; j <= ncom; j++)
{
xt[j] = pcom[j] + x * xicom[j];
}
return func2(xt, ncom);
}
int frprmn_iter;
double frprmn_fret;
void frprmn(double p[], int n, double ftol, int iter, double fret)
{
int j, its, itmax = 200;
double eps = 0.0000000001;
double g[] = new double[551];
double h[] = new double[51];
double xi[] = new double[51];
double fp = func2(p, n);
dfunc(p, xi);
for (j = 1; j <= n; j++)
{
g[j] = -xi[j];
h[j] = g[j];
xi[j] = h[j];
}
for (its = 1; its <= itmax; its++)
{
iter = its;
linmin(p, xi, n, fret);
fret = linmin_fret;
if (2.0 * Math.abs(fret - fp) <= ftol * (Math.abs(fret) + Math.abs(fp) + eps))
{
break;
}
fp = func2(p, n);
dfunc(p, xi);
double gg = 0.0;
double dgg = 0.0;
for (j = 1; j <= n; j++)
{
gg = gg + g[j] * g[j];
//dgg = dgg + xi(j) * xi(j) //polak-ribiere 法
dgg = dgg + (xi[j] + g[j]) * xi[j]; //fletcher-reeves 法
}
if (gg == 0.0)
{
break;
}
double gam = dgg / gg;
for (j = 1; j <= n; j++)
{
g[j] = -xi[j];
h[j] = g[j] + gam * h[j];
xi[j] = h[j];
}
}
if (iter > itmax)
{
System.out.println("frpr maximum iterations exceeded");
System.exit(1);
}
frprmn_iter = iter;
frprmn_fret = fret;
}
double bessj0(double x)
{
double p1, p2, p3, p4, p5, q1, q2, q3, q4, q5;
double r1, r2, r3, r4, r5, r6, s1, s2, s3, s4, s5, s6;
double y, bbb, ccc, aaa, temp, eee, ddd, ax, xx, z;
p1 = 1.0; p2 = -0.001098628627;
p3 = 0.00002734510407; p4 = -0.000002073370639;
p5 = 2.093887211e-07;
q1 = -0.01562499995; q2 = 0.0001430488765;
q3 = -0.000006911147651; q4 = 7.621095161e-07;
q5 = -9.34945152e-08;
r1 = 57568490574.0; r2 = -.0;
r3 = 651619640.7; r4 = -11214424.18;
r5 = 77392.33017; r6 = -184.9052456;
s1 = 57568490411.0; s2 = 1029532985.0;
s3 = 9494680.718; s4 = 59272.64853;
s5 = 267.8532712; s6 = 1.0;
if (Math.abs(x) < 8.0)
{
y = x * x;
bbb = y * (r4 + y * (r5 + y * r6));
aaa = r1 + y * (r2 + y * (r3 + bbb));
ccc = y * (s3 + y * (s4 + y * (s5 + y * s6)));
temp = aaa / (s1 + y * (s2 + ccc));
}
else
{
ax = Math.abs(x);
z = 8.0 / ax;
y = z * z;
xx = ax - 0.785398164;
ccc = y * (p3 + y * (p4 + y * p5));
aaa = p1 + y * (p2 + ccc);
ddd = y * (q3 + y * (q4 + y * q5));
eee = z * Math.sin(xx) * (q1 + y * (q2 + ddd));
temp = Math.sqrt(0.636619772 / ax) * (Math.cos(xx) * aaa - eee);
}
return temp;
}
double bessj1(double x)
{
double p1, p2, p3, p4, p5, q1, q2, q3, q4, q5;
double r1, r2, r3, r4, r5, r6, s1, s2, s3, s4, s5, s6;
double bbb, ccc, aaa, temp, ax, xx, z, y;
r1 = 72362614232.0; r2 = -7895059235.0;
r3 = 242396853.1; r4 = -2972611.439;
r5 = 15704.4826; r6 = -30.16036606;
s1 = 144725228442.0; s2 = 2300535178.0;
s3 = 18583304.74; s4 = 99447.43394;
s5 = 376.9991397; s6 = 1.0;
p1 = 1.0; p2 = 0.00183105;
p3 = -0.00003516396496; p4 = 0.000002457520174;
p5 = -0.000000240337019;
q1 = 0.04687499995; q2 = -0.0002002690873;
q3 = 0.000008449199096; q4 = -0.00000088228987;
q5 = 0.000000105787412;
if (Math.abs(x) < 8.0 )
{
y = x * x;
aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
bbb = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
temp = x * aaa / bbb;
}
else
{
ax = Math.abs(x);
z = 8.0 / ax;
y = z * z;
xx = ax - 2.356194491;
aaa = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
bbb = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
ccc = Math.sqrt(0.636619772 / ax);
temp = ccc * (Math.cos(xx) * aaa - z * Math.sin(xx) * bbb * sgn(x));
}
return temp;
}
int sgn(double x)
{
int temp = 0;
if (x > 0) temp = 1;
if (x == 0) temp = 0;
if (x < 0) temp = -1;
return temp;
}
int ncom;
double pcom[] = new double[51];
double xicom[] = new double[51];
double linmin_fret;
void linmin(double p[], double xi[], int n, double fret)
{
int j;
double tol = 0.0001;
ncom = n;
for (j = 1; j <= n; j++)
{
pcom[j] = p[j];
xicom[j] = xi[j];
}
double fa, fx, fb, bx, ax = 0.0;
fa = 0; fb = 0; fx = 0; bx = 0;
double xmin, xx = 1.0;
xmin = 0;
mnbrak(ax, xx, bx, fa, fx, fb);
ax = mnbrak_ax; xx = mnbrak_bx; bx = mnbrak_cx; fa = mnbrak_fa; fx = mnbrak_fb; fb = mnbrak_fc;
fret = brent(ax, xx, bx, tol, xmin);
xmin = brent_xmin;
for (j = 1; j <= n; j++)
{
xi[j] = xmin * xi[j];
p[j] = p[j] + xi[j];
}
linmin_fret = fret;
}
double brent_xmin;
double brent(double ax, double bx, double cx, double tol, double xmin)
{
boolean done;
int iter,itmax = 100;
double d, fu, r, q, p, xm, tol1, tol2, a, b, cgold = 0.381966;
double u, etemp, dum, v, w, x, e, fx, fv1, fw, zeps = 0.0000000001;
d = 0;
a = ax;
if (cx < ax)
{
a = cx;
}
b = ax;
if (cx > ax)
{
b = cx;
}
v = bx;
w = v;
x = v;
e = 0.0;
fx = func(x);
fv1 = fx;
fw = fx;
for (iter = 1; iter <= itmax; iter++)
{
xm = 0.5 * (a + b);
tol1 = tol * Math.abs(x) + zeps;
tol2 = 2.0 * tol1;
if (Math.abs(x - xm) <= tol2 - 0.5 * (b - a))
{
break;
}
done = true;
if (Math.abs(e) > tol1)
{
r = (x - w) * (fx - fv1);
q = (x - v) * (fx - fw);
p = (x - v) * q - (x - w) * r;
q = 2.0 * (q - r);
if (q > 0.0)
{
p = -p;
}
q = Math.abs(q);
etemp = e;
e = d;
dum = Math.abs(0.5 * q * etemp);
if (Math.abs(p) < dum && p > q * (a - x) && p < q * (b - x))
{
d = p / q;
u = x + d;
if (u - a < tol2 || b - u < tol2)
{
d = Math.abs(tol1) * sgn(xm - x);
}
done = false;
}
}
if (done)
{
if (x >= xm)
{
e = a - x;
}
else
{
e = b - x;
}
d = cgold * e;
}
if (Math.abs(d) >= tol1)
{
u = x + d;
}
else
{
u = x + Math.abs(tol1) * sgn(d);
}
fu = func(u);
if (fu <= fx)
{
if (u >= x)
{
a = x;
}
上一篇:
d11r8.java
下一篇:
那些令你为之触动的好句子,感慨万分