📄 d14r4f.java
字号:
public class d14r4F
{
int kount; int kmax; double dxsav;
int odeint_nvar;
double odeint_x1, odeint_x2, odeint_eps, odeint_h1, odeint_hmin, odeint_nok, odeint_nbad;
void odeint(double ystart[], int nvar, double x1, double x2, double eps, double h1, double hmin, double nok, double nbad, double xp[], double yp[][])
{
int i, nstp, maxstp;
double yscal[] = new double[11];
double y[] = new double[11];
double dydx[] = new double[11];
double two, zero, tiny, x, h, xsav, hdid, hnext;
xsav = 0;
maxstp = 10000;
two = 2.0;
zero = 0.0;
tiny = 1e-30;
x = x1;
h = Math.abs(h1) * sgn(x2 - x1);
nok = 0;
nbad = 0;
kount = 0;
hdid = 0.0;
hnext = 0.0;
for (i = 1; i <= nvar; i++)
{
y[i] = ystart[i];
}
if (kmax > 0)
xsav = x - dxsav * two;
for (nstp = 1; nstp <= maxstp; nstp++)
{
derivs(x, y, dydx);
for (i = 1; i <= nvar; i++)
{
yscal[i] = Math.abs(y[i]) + Math.abs(h * dydx[i]) + tiny;
}
if (kmax > 0 )
{
if (Math.abs(x - xsav) > Math.abs(dxsav))
{
if (kount < kmax - 1)
{
kount = kount + 1;
xp[kount] = x;
for (i = 1; i <= nvar; i++)
{
yp[i][kount] = y[i];
}
xsav = x;
}
}
}
if( (x + h - x2) * (x + h - x1) > zero )
{
h = x2 - x;
}
rkqc(y, dydx, nvar, x, h, eps, yscal, hdid, hnext);
x = rkqc_x; hdid = rkqc_hdid; hnext = rkqc_hnext;
if (hdid == h)
{
nok = nok + 1;
}
else
{
nbad = nbad + 1;
}
if ((x - x2) * (x2 - x1) >= zero )
{
for (i = 1; i <= nvar; i++)
{
ystart[i] = y[i];
}
if (kmax != 0)
{
kount = kount + 1;
xp[kount] = x;
for (i = 1; i <= nvar; i++)
{
yp[i][kount] = y[i];
}
}
for (i = 1; i <= 10; i++)
{
dydx[i] = 0.0;
y[i] = 0.0;
yscal[i] = 0.0;
}
odeint_x1 = x1; odeint_x2 = x2; odeint_eps = eps; odeint_h1 = h1; odeint_hmin = hmin;
odeint_nok = nok; odeint_nbad = nbad; odeint_nvar = nvar;
return;
}
if (Math.abs(hnext) < hmin)
{
System.out.println(" stepsize smaller than minimum.");
System.exit(1);
}
h = hnext;
}
System.out.println(" too many steps.");
System.exit(1);
}
double rkqc_x, rkqc_hdid, rkqc_hnext;
void rkqc(double y[], double dydx[], int n, double x, double htry, double eps, double yscal[], double hdid, double hnext)
{
int i, flag;
double one, safety, errcon, fcor;
double ytemp[] = new double[11];
double ysav[] = new double[11];
double dysav[] = new double[11];
double xsav, h, hh, pgrow, pshrnk, errmax;
one = 1.0;
safety = 0.9;
errcon = 0.0006;
fcor = 0.066667;
pgrow = -0.2;
pshrnk = -0.25;
xsav = x;
for (i = 1; i <= n; i++)
{
ysav[i] = y[i];
dysav[i] = dydx[i];
}
h = htry;
do
{
hh = 0.5 * h;
rk4(ysav, dysav, n, xsav, hh, ytemp);
x = xsav + hh;
derivs(x, ytemp, dydx);
rk4(ytemp, dydx, n, x, hh, y);
x = xsav + h;
if( x == xsav)
{
System.out.println(" stepsize not significant in rkqc.");
System.exit(1);
}
rk4(ysav, dysav, n, xsav, h, ytemp);
errmax = 0.0;
for (i = 1; i <= n; i++)
{
ytemp[i] = y[i] - ytemp[i];
if (errmax < Math.abs(ytemp[i] / yscal[i]))
{
errmax = Math.abs(ytemp[i] / yscal[i]);
}
}
errmax = errmax / eps;
if (errmax > one)
{
h = safety * h * Math.pow(errmax, pshrnk);
flag = 1;
}
else
{
hdid = h;
if( errmax > errcon )
{
hnext = safety * h * Math.pow(errmax , pgrow);
}
else
{
hnext = 4.0 * h;
}
flag = 0;
}
}
while( flag == 1);
for (i = 1; i <= n; i++)
{
y[i] = y[i] + ytemp[i] * fcor;
}
for (i = 1; i <= 10; i++)
{
dysav[i] = 0.0;
ysav[i] = 0.0;
ytemp[i] = 0.0;
}
rkqc_x = x; rkqc_hdid = hdid; rkqc_hnext = hnext;
}
void rk4(double y[], double dydx[], int n, double x, double h, double yout[])
{
int i;
double yt[] = new double[11];
double dyt[] = new double[11];
double dym[] = new double[11];
double hh, h6, xh;
hh = h * 0.5;
h6 = h / 6.0;
xh = x + hh;
for( i = 1; i <= n; i++)
{
yt[i] = y[i] + hh * dydx[i];
}
derivs(xh, yt, dyt);
for( i = 1; i <= n; i++)
{
yt[i] = y[i] + hh * dyt[i];
}
derivs(xh, yt, dym);
for( i = 1; i <= n; i++)
{
yt[i] = y[i] + h * dym[i];
dym[i] = dyt[i] + dym[i];
}
derivs(x + h, yt, dyt);
for (i = 1; i <= n; i++)
{
yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
}
for(i=1; i <= 10; i++)
{
dym[i]=0.0;
dyt[i]=0.0;
yt[i]=0.0;
}
}
void derivs(double x,double y1[],double dydx[])
{
dydx[1] = -y1[2];
dydx[2] = y1[1] - (1.0 / x) * y1[2];
dydx[3] = y1[2] - (2.0 / x) * y1[3];
dydx[4] = y1[3] - (3.0 / x) * y1[4];
}
double bessj(int n, double x)
{
int iacc, j, m;
double temp, bigno, bigni, ax, tox, bj, bjp, bjm, besj, jsum, sum;
iacc = 40;
bigno = 10000000000.0;
bigni = 0.0000000001;
if (n < 2 )
{
System.out.println("bad argument n in bassj");
System.exit(1);
}
ax = Math.abs(x);
if (ax == 0)
{
temp = 0.0;
}
else if (ax > (float)n)
{
tox = 2.0 / ax;
bjm = bessj0(ax);
bj = bessj1(ax);
for (j = 1; j <= n - 1; j++)
{
bjp = j * tox * bj - bjm;
bjm = bj;
bj = bjp;
}
temp = bj;
}
else
{
tox = 2.0 / ax;
m = 2 * (int)(((n + (int)(Math.sqrt(iacc * n)))) / 2);
besj = 0.0;
jsum = 0;
sum = 0.0;
bjp = 0.0;
bj = 1.0;
for (j = m; j >= 1; j--)
{
bjm = j * tox * bj - bjp;
bjp = bj;
bj = bjm;
if (Math.abs(bj) > bigno)
{
bj = bj * bigni;
bjp = bjp * bigni;
besj = besj * bigni;
sum = sum * bigni;
}
if (jsum != 0) sum = sum + bj;
jsum = 1 - jsum;
if (j == n) besj = bjp;
}
sum = 2.0 * sum - bj;
temp = besj / sum;
}
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;
}
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 = -13362590354.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;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -