16 #define ODE_NVAR_MAX 256
17 static int RK4Solve (
double *,
int,
double,
double,
18 void (*rhs)(
double,
double *,
double *));
19 static int CK45Solve (
double *,
int,
double,
double *,
20 void (*rhs)(
double,
double *,
double *),
double);
23 void ODE_Solve(
double *y0,
int nvar,
double xbeg,
double xend,
24 double dx,
void (*rhs)(
double,
double *,
double *),
int method)
53 nsteps = (int)((xend - xbeg)/dx);
79 err =
CK45Solve (y0, nvar, x0, &dx, rhs, 1.e-9);
84 if (fabs(dx/dx0) < 1.e-4) {
85 print (
"! ODE_Solve: step size too small \n");
89 print (
"! ODE_Solve: too many steps\n");
90 print (
"! xbeg = %f, xend = %f\n",xbeg, xend);
96 print1 (
"! ODE_Solve: integration method unknown\n");
104 double dx,
void (*rhs)(
double,
double *,
double *))
118 for (nv = 0; nv < nvar; nv++) y1[nv] = y0[nv] + 0.5*dx*k1[nv];
122 rhs (x0 + 0.5*dx, y1, k2);
123 for (nv = 0; nv < nvar; nv++) y1[nv] = y0[nv] + 0.5*dx*k2[nv];
127 rhs (x0 + 0.5*dx, y1, k3);
128 for (nv = 0; nv < nvar; nv++) y1[nv] = y0[nv] + dx*k3[nv];
132 rhs (x0 + dx, y1, k4);
133 for (nv = 0; nv < nvar; nv++)
134 y0[nv] += dx*(k1[nv] + 2.0*(k2[nv] + k3[nv]) + k4[nv])/6.0;
141 double *dx0,
void (*rhs)(
double,
double *,
double *),
double tol)
148 int i,
n, rhs_err, isgn, useZ=0;
150 double scrh, err, dx_shrink, dx_grow;
153 const double a2 = 0.2;
154 const double a3 = 0.3;
155 const double a4 = 0.6;
156 const double a5 = 1.0;
157 const double a6 = 0.875;
159 const double c1 = 37.0/378.0;
160 const double c2 = 0.0;
161 const double c3 = 250.0/621.0;
162 const double c4 = 125.0/594.0;
163 const double c5 = 0.0;
164 const double c6 = 512.0/1771.0;
166 const double cs1 = 2825.0/27648.0;
167 const double cs2 = 0.0;
168 const double cs3 = 18575.0/48384.0;
169 const double cs4 = 13525.0/55296.0;
170 const double cs5 = 277.0/14336.0;
171 const double cs6 = 0.25;
173 const double b21 = 0.2;
174 const double b31 = 3.0/40.0 , b32 = 9.0/40.0;
175 const double b41 = 0.3 , b42 = -0.9, b43 = 1.2;
176 const double b51 = -11.0/54.0, b52 = 2.5, b53 = -70.0/27.0, b54 = 35.0/27.0;
177 const double b61 = 1631.0/55296.0, b62 = 175.0/512.0, b63 = 575.0/13824.0, b64 = 44275.0/110592.0, b65 = 253.0/4096.0;
188 if (dx < 0.0) isgn = -1;
193 for (n = 0; n < nvar; n++) y1[n] = y0[n];
197 for (n = 0; n < nvar; n++) ys[n] = y1[n] + dx*b21*k1[n];
199 rhs (x0 + a2*dx, ys, k2);
200 for (n = 0; n < nvar; n++) ys[n] = y1[n] + dx*(b31*k1[n] + b32*k2[n]);
202 rhs (x0 + a3*dx, ys, k3);
203 for (n = 0; n < nvar; n++) {
204 ys[
n] = y1[
n] + dx*(b41*k1[
n] + b42*k2[
n] + b43*k3[
n]);
207 rhs (x0 + a4*dx, ys, k4);
208 for (n = 0; n < nvar; n++) {
209 ys[
n] = y1[
n] + dx*(b51*k1[
n] + b52*k2[
n] + b53*k3[
n] + b54*k4[
n]);
212 rhs (x0 + a5*dx, ys, k5);
213 for (n = 0; n < nvar; n++) {
214 ys[
n] = y1[
n] + dx*(b61*k1[
n] + b62*k2[
n] + b63*k3[
n]
215 + b64*k4[
n] + b65*k5[
n]);
218 rhs (x0 + a6*dx, ys, k6);
222 for (n = 0; n < nvar; n++) {
223 y5th[
n] = y1[
n] + dx*(c1*k1[
n] + c2*k2[
n] + c3*k3[
n]
224 + c4*k4[
n] + c5*k5[
n] + c6*k6[
n]);
225 y4th[
n] = y1[
n] + dx*(cs1*k1[
n] + cs2*k2[
n] + cs3*k3[
n]
226 + cs4*k4[
n] + cs5*k5[
n] + cs6*k6[
n]);
238 for (n = 0; n < nvar; n++) {
239 scrh = fabs(y5th[n] - y4th[n])/(1.0 + fabs(y1[n]));
240 err =
MAX(err, scrh);
247 err =
MAX(err, 1.e-18);
251 dx_grow = 0.9*fabs(dx)*pow(err, -0.2);
252 dx_grow =
MIN(dx_grow, 5.0*fabs(dx));
255 for (n = 0; n < nvar; n++) y0[n] = y1[n] = y5th[n];
260 dx_shrink = 0.9*fabs(dx)*pow(err, -0.25);
261 dx =
MAX(dx_shrink, 0.05*fabs(dx));
void print1(const char *fmt,...)
void ODE_Solve(double *y0, int nvar, double xbeg, double xend, double dx, void(*rhs)(double, double *, double *), int method)
static int RK4Solve(double *, int, double, double, void(*rhs)(double, double *, double *))
void print(const char *fmt,...)
static int CK45Solve(double *, int, double, double *, void(*rhs)(double, double *, double *), double)
#define QUIT_PLUTO(e_code)