Bonsoir,
J'ai un vrai souci avec numerical recipes 3.
C'est comme si ce code ne respectait pas les bonnes pratiques cpp !
Si je fais un
#include "stepperdopr853.h"
dans 2 fihiers je recupere :
multiple definition of `Dopr853_constants::
Voila ce que j'ai trouve
...no other file contains
#include "stepperdopr853.h" Do you have a header file that includes stepperdopr853.h? Is that header included from more than one .cc file? (For example: Does bratoket.h include stepperdopr853.h and is bratoket.h included in more than one of your .cc files?)
#ifndef ... doesn't helpThat won't protect against including the same header file from other .cc files.
Each .cc file is compiled separately and a symbol table is created for that object file. When the linker puts all of the object files together it merges the symbol tables. If the same symbol is defined in more than one of the object file tables, the linker reports an error and no executable can be created.
Bottom line: With files from the Numerical Recipes Version 3 distribution you can't include header files like stepperdopr853.h from more than one .cc file since the header file actually declares memory storage.
It "just happens" that, since stepperdopr5.h just consists of some templates and doesn't contain any data definitions that would result in problems like this, you actually can get away with including it from more than one .cc file. But that's not true for many other NR3 header files.
If you are wondering why the package was put together that way:
See responses from one of the authors here:
http://www.nr.com/forum/showthread.php?t=1332 and here:
http://www.nr.com/forum/showthread.php?t=1260
Helas, ces forums n'existent plus
J'attache le fichier qui pose problem.
Quelle solution je peux apporter ?
faire 2 fichiers
L'un avec les constants ?
Je suppose que ces les constantes static qui font le probleme.
Merci d'avance
===========================================================================
Patrick DUPRÉ | | email: pdupre@???
===========================================================================
#pragma once
#include <nr3.h>
#include "stepper.h"
struct Dopr853_constants {
static const Doub c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c14, c15, c16,
b1, b6, b7, b8, b9, b10, b11, b12, bhh1, bhh2, bhh3,
er1, er6, er7, er8, er9, er10, er11, er12,
a21, a31, a32, a41, a43, a51, a53, a54, a61, a64, a65, a71, a74, a75, a76,
a81, a84, a85, a86, a87, a91, a94, a95, a96, a97, a98, a101, a104, a105,
a106, a107, a108, a109, a111, a114, a115, a116, a117, a118, a119, a1110,
a121, a124, a125, a126, a127, a128, a129, a1210, a1211, a141, a147, a148,
a149, a1410, a1411, a1412, a1413, a151, a156, a157, a158, a1511, a1512,
a1513, a1514, a161, a166, a167, a168, a169, a1613, a1614, a1615,
d41, d46, d47, d48, d49, d410, d411, d412, d413, d414, d415, d416, d51, d56,
d57, d58, d59, d510, d511, d512, d513, d514, d515, d516, d61, d66, d67, d68,
d69, d610, d611, d612, d613, d614, d615, d616, d71, d76, d77, d78, d79,
d710, d711, d712, d713, d714, d715, d716;
};
const Doub Dopr853_constants::c2 = 0.526001519587677318785587544488e-01;
const Doub Dopr853_constants::c3 = 0.789002279381515978178381316732e-01;
const Doub Dopr853_constants::c4 = 0.118350341907227396726757197510e+00;
const Doub Dopr853_constants::c5 = 0.281649658092772603273242802490e+00;
const Doub Dopr853_constants::c6 = 0.333333333333333333333333333333e+00;
const Doub Dopr853_constants::c7 = 0.25e+00;
const Doub Dopr853_constants::c8 = 0.307692307692307692307692307692e+00;
const Doub Dopr853_constants::c9 = 0.651282051282051282051282051282e+00;
const Doub Dopr853_constants::c10 = 0.6e+00;
const Doub Dopr853_constants::c11 = 0.857142857142857142857142857142e+00;
const Doub Dopr853_constants::c14 = 0.1e+00;
const Doub Dopr853_constants::c15 = 0.2e+00;
const Doub Dopr853_constants::c16 = 0.777777777777777777777777777778e+00;
const Doub Dopr853_constants::b1 = 5.42937341165687622380535766363e-2;
const Doub Dopr853_constants::b6 = 4.45031289275240888144113950566e0;
const Doub Dopr853_constants::b7 = 1.89151789931450038304281599044e0;
const Doub Dopr853_constants::b8 = -5.8012039600105847814672114227e0;
const Doub Dopr853_constants::b9 = 3.1116436695781989440891606237e-1;
const Doub Dopr853_constants::b10 = -1.52160949662516078556178806805e-1;
const Doub Dopr853_constants::b11 = 2.01365400804030348374776537501e-1;
const Doub Dopr853_constants::b12 = 4.47106157277725905176885569043e-2;
const Doub Dopr853_constants::bhh1 = 0.244094488188976377952755905512e+00;
const Doub Dopr853_constants::bhh2 = 0.733846688281611857341361741547e+00;
const Doub Dopr853_constants::bhh3 = 0.220588235294117647058823529412e-01;
const Doub Dopr853_constants::er1 = 0.1312004499419488073250102996e-01;
const Doub Dopr853_constants::er6 = -0.1225156446376204440720569753e+01;
const Doub Dopr853_constants::er7 = -0.4957589496572501915214079952e+00;
const Doub Dopr853_constants::er8 = 0.1664377182454986536961530415e+01;
const Doub Dopr853_constants::er9 = -0.3503288487499736816886487290e+00;
const Doub Dopr853_constants::er10 = 0.3341791187130174790297318841e+00;
const Doub Dopr853_constants::er11 = 0.8192320648511571246570742613e-01;
const Doub Dopr853_constants::er12 = -0.2235530786388629525884427845e-01;
const Doub Dopr853_constants::a21 = 5.26001519587677318785587544488e-2;
const Doub Dopr853_constants::a31 = 1.97250569845378994544595329183e-2;
const Doub Dopr853_constants::a32 = 5.91751709536136983633785987549e-2;
const Doub Dopr853_constants::a41 = 2.95875854768068491816892993775e-2;
const Doub Dopr853_constants::a43 = 8.87627564304205475450678981324e-2;
const Doub Dopr853_constants::a51 = 2.41365134159266685502369798665e-1;
const Doub Dopr853_constants::a53 = -8.84549479328286085344864962717e-1;
const Doub Dopr853_constants::a54 = 9.24834003261792003115737966543e-1;
const Doub Dopr853_constants::a61 = 3.7037037037037037037037037037e-2;
const Doub Dopr853_constants::a64 = 1.70828608729473871279604482173e-1;
const Doub Dopr853_constants::a65 = 1.25467687566822425016691814123e-1;
const Doub Dopr853_constants::a71 = 3.7109375e-2;
const Doub Dopr853_constants::a74 = 1.70252211019544039314978060272e-1;
const Doub Dopr853_constants::a75 = 6.02165389804559606850219397283e-2;
const Doub Dopr853_constants::a76 = -1.7578125e-2;
const Doub Dopr853_constants::a81 = 3.70920001185047927108779319836e-2;
const Doub Dopr853_constants::a84 = 1.70383925712239993810214054705e-1;
const Doub Dopr853_constants::a85 = 1.07262030446373284651809199168e-1;
const Doub Dopr853_constants::a86 = -1.53194377486244017527936158236e-2;
const Doub Dopr853_constants::a87 = 8.27378916381402288758473766002e-3;
const Doub Dopr853_constants::a91 = 6.24110958716075717114429577812e-1;
const Doub Dopr853_constants::a94 = -3.36089262944694129406857109825e0;
const Doub Dopr853_constants::a95 = -8.68219346841726006818189891453e-1;
const Doub Dopr853_constants::a96 = 2.75920996994467083049415600797e1;
const Doub Dopr853_constants::a97 = 2.01540675504778934086186788979e1;
const Doub Dopr853_constants::a98 = -4.34898841810699588477366255144e1;
const Doub Dopr853_constants::a101 = 4.77662536438264365890433908527e-1;
const Doub Dopr853_constants::a104 = -2.48811461997166764192642586468e0;
const Doub Dopr853_constants::a105 = -5.90290826836842996371446475743e-1;
const Doub Dopr853_constants::a106 = 2.12300514481811942347288949897e1;
const Doub Dopr853_constants::a107 = 1.52792336328824235832596922938e1;
const Doub Dopr853_constants::a108 = -3.32882109689848629194453265587e1;
const Doub Dopr853_constants::a109 = -2.03312017085086261358222928593e-2;
const Doub Dopr853_constants::a111 = -9.3714243008598732571704021658e-1;
const Doub Dopr853_constants::a114 = 5.18637242884406370830023853209e0;
const Doub Dopr853_constants::a115 = 1.09143734899672957818500254654e0;
const Doub Dopr853_constants::a116 = -8.14978701074692612513997267357e0;
const Doub Dopr853_constants::a117 = -1.85200656599969598641566180701e1;
const Doub Dopr853_constants::a118 = 2.27394870993505042818970056734e1;
const Doub Dopr853_constants::a119 = 2.49360555267965238987089396762e0;
const Doub Dopr853_constants::a1110 = -3.0467644718982195003823669022e0;
const Doub Dopr853_constants::a121 = 2.27331014751653820792359768449e0;
const Doub Dopr853_constants::a124 = -1.05344954667372501984066689879e1;
const Doub Dopr853_constants::a125 = -2.00087205822486249909675718444e0;
const Doub Dopr853_constants::a126 = -1.79589318631187989172765950534e1;
const Doub Dopr853_constants::a127 = 2.79488845294199600508499808837e1;
const Doub Dopr853_constants::a128 = -2.85899827713502369474065508674e0;
const Doub Dopr853_constants::a129 = -8.87285693353062954433549289258e0;
const Doub Dopr853_constants::a1210 = 1.23605671757943030647266201528e1;
const Doub Dopr853_constants::a1211 = 6.43392746015763530355970484046e-1;
const Doub Dopr853_constants::a141 = 5.61675022830479523392909219681e-2;
const Doub Dopr853_constants::a147 = 2.53500210216624811088794765333e-1;
const Doub Dopr853_constants::a148 = -2.46239037470802489917441475441e-1;
const Doub Dopr853_constants::a149 = -1.24191423263816360469010140626e-1;
const Doub Dopr853_constants::a1410 = 1.5329179827876569731206322685e-1;
const Doub Dopr853_constants::a1411 = 8.20105229563468988491666602057e-3;
const Doub Dopr853_constants::a1412 = 7.56789766054569976138603589584e-3;
const Doub Dopr853_constants::a1413 = -8.298e-3;
const Doub Dopr853_constants::a151 = 3.18346481635021405060768473261e-2;
const Doub Dopr853_constants::a156 = 2.83009096723667755288322961402e-2;
const Doub Dopr853_constants::a157 = 5.35419883074385676223797384372e-2;
const Doub Dopr853_constants::a158 = -5.49237485713909884646569340306e-2;
const Doub Dopr853_constants::a1511 = -1.08347328697249322858509316994e-4;
const Doub Dopr853_constants::a1512 = 3.82571090835658412954920192323e-4;
const Doub Dopr853_constants::a1513 = -3.40465008687404560802977114492e-4;
const Doub Dopr853_constants::a1514 = 1.41312443674632500278074618366e-1;
const Doub Dopr853_constants::a161 = -4.28896301583791923408573538692e-1;
const Doub Dopr853_constants::a166 = -4.69762141536116384314449447206e0;
const Doub Dopr853_constants::a167 = 7.68342119606259904184240953878e0;
const Doub Dopr853_constants::a168 = 4.06898981839711007970213554331e0;
const Doub Dopr853_constants::a169 = 3.56727187455281109270669543021e-1;
const Doub Dopr853_constants::a1613 = -1.39902416515901462129418009734e-3;
const Doub Dopr853_constants::a1614 = 2.9475147891527723389556272149e0;
const Doub Dopr853_constants::a1615 = -9.15095847217987001081870187138e0;
const Doub Dopr853_constants::d41 = -0.84289382761090128651353491142e+01;
const Doub Dopr853_constants::d46 = 0.56671495351937776962531783590e+00;
const Doub Dopr853_constants::d47 = -0.30689499459498916912797304727e+01;
const Doub Dopr853_constants::d48 = 0.23846676565120698287728149680e+01;
const Doub Dopr853_constants::d49 = 0.21170345824450282767155149946e+01;
const Doub Dopr853_constants::d410 = -0.87139158377797299206789907490e+00;
const Doub Dopr853_constants::d411 = 0.22404374302607882758541771650e+01;
const Doub Dopr853_constants::d412 = 0.63157877876946881815570249290e+00;
const Doub Dopr853_constants::d413 = -0.88990336451333310820698117400e-01;
const Doub Dopr853_constants::d414 = 0.18148505520854727256656404962e+02;
const Doub Dopr853_constants::d415 = -0.91946323924783554000451984436e+01;
const Doub Dopr853_constants::d416 = -0.44360363875948939664310572000e+01;
const Doub Dopr853_constants::d51 = 0.10427508642579134603413151009e+02;
const Doub Dopr853_constants::d56 = 0.24228349177525818288430175319e+03;
const Doub Dopr853_constants::d57 = 0.16520045171727028198505394887e+03;
const Doub Dopr853_constants::d58 = -0.37454675472269020279518312152e+03;
const Doub Dopr853_constants::d59 = -0.22113666853125306036270938578e+02;
const Doub Dopr853_constants::d510 = 0.77334326684722638389603898808e+01;
const Doub Dopr853_constants::d511 = -0.30674084731089398182061213626e+02;
const Doub Dopr853_constants::d512 = -0.93321305264302278729567221706e+01;
const Doub Dopr853_constants::d513 = 0.15697238121770843886131091075e+02;
const Doub Dopr853_constants::d514 = -0.31139403219565177677282850411e+02;
const Doub Dopr853_constants::d515 = -0.93529243588444783865713862664e+01;
const Doub Dopr853_constants::d516 = 0.35816841486394083752465898540e+02;
const Doub Dopr853_constants::d61 = 0.19985053242002433820987653617e+02;
const Doub Dopr853_constants::d66 = -0.38703730874935176555105901742e+03;
const Doub Dopr853_constants::d67 = -0.18917813819516756882830838328e+03;
const Doub Dopr853_constants::d68 = 0.52780815920542364900561016686e+03;
const Doub Dopr853_constants::d69 = -0.11573902539959630126141871134e+02;
const Doub Dopr853_constants::d610 = 0.68812326946963000169666922661e+01;
const Doub Dopr853_constants::d611 = -0.10006050966910838403183860980e+01;
const Doub Dopr853_constants::d612 = 0.77771377980534432092869265740e+00;
const Doub Dopr853_constants::d613 = -0.27782057523535084065932004339e+01;
const Doub Dopr853_constants::d614 = -0.60196695231264120758267380846e+02;
const Doub Dopr853_constants::d615 = 0.84320405506677161018159903784e+02;
const Doub Dopr853_constants::d616 = 0.11992291136182789328035130030e+02;
const Doub Dopr853_constants::d71 = -0.25693933462703749003312586129e+02;
const Doub Dopr853_constants::d76 = -0.15418974869023643374053993627e+03;
const Doub Dopr853_constants::d77 = -0.23152937917604549567536039109e+03;
const Doub Dopr853_constants::d78 = 0.35763911791061412378285349910e+03;
const Doub Dopr853_constants::d79 = 0.93405324183624310003907691704e+02;
const Doub Dopr853_constants::d710 = -0.37458323136451633156875139351e+02;
const Doub Dopr853_constants::d711 = 0.10409964950896230045147246184e+03;
const Doub Dopr853_constants::d712 = 0.29840293426660503123344363579e+02;
const Doub Dopr853_constants::d713 = -0.43533456590011143754432175058e+02;
const Doub Dopr853_constants::d714 = 0.96324553959188282948394950600e+02;
const Doub Dopr853_constants::d715 = -0.39177261675615439165231486172e+02;
const Doub Dopr853_constants::d716 = -0.14972683625798562581422125276e+03;
template<class D>
struct StepperDopr853 : StepperBase, Dopr853_constants {
typedef D Dtype;
VecDoub yerr2;
VecDoub k2, k3, k4, k5, k6, k7, k8, k9, k10;
VecDoub rcont1, rcont2, rcont3, rcont4, rcont5, rcont6, rcont7, rcont8;
StepperDopr853(VecDoub_IO &yy, VecDoub_IO &dydxx, Doub &xx,
const Doub atoll, const Doub rtoll, bool dens);
void step(const Doub htry, D &derivs);
void dy(const Doub h, D &derivs);
void prepare_dense(const Doub h, VecDoub_I &dydxnew, D &derivs);
Doub dense_out(const Int i, const Doub x, const Doub h);
Doub error(const Doub h);
struct Controller {
Doub hnext, errold;
bool reject;
Controller();
bool success(const Doub err, Doub &h);
};
Controller con;
};
template<class D>
StepperDopr853<D>::StepperDopr853(VecDoub_IO &yy, VecDoub_IO &dydxx, Doub &xx,
const Doub atoll, const Doub rtoll, bool dens) :
StepperBase(yy, dydxx, xx, atoll, rtoll, dens), yerr2(n), k2(n), k3(n), k4(n),
k5(n), k6(n), k7(n), k8(n), k9(n), k10(n), rcont1(n), rcont2(n), rcont3(n),
rcont4(n), rcont5(n), rcont6(n), rcont7(n), rcont8(n) {
EPS = std::numeric_limits<Doub>::epsilon();
}
template<class D>
void StepperDopr853<D>::step(const Doub htry, D &derivs) {
VecDoub dydxnew(n);
Doub h = htry;
for (;;) {
dy(h, derivs);
Doub err = error(h);
if (con.success(err, h)) break;
if (abs(h) <= abs(x) * EPS) throw("stepsize underflow in StepperDopr853");
}
derivs(x + h, yout, dydxnew);
if (dense)
prepare_dense(h, dydxnew, derivs);
dydx = dydxnew;
y = yout;
xold = x;
x += (hdid = h);
hnext = con.hnext;
}
template<class D>
void StepperDopr853<D>::dy(const Doub h, D &derivs) {
VecDoub ytemp(n);
Int i;
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * a21 * dydx[i];
derivs(x + c2 * h, ytemp, k2);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a31 * dydx[i] + a32 * k2[i]);
derivs(x + c3 * h, ytemp, k3);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a41 * dydx[i] + a43 * k3[i]);
derivs(x + c4 * h, ytemp, k4);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a51 * dydx[i] + a53 * k3[i] + a54 * k4[i]);
derivs(x + c5 * h, ytemp, k5);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a61 * dydx[i] + a64 * k4[i] + a65 * k5[i]);
derivs(x + c6 * h, ytemp, k6);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a71 * dydx[i] + a74 * k4[i] + a75 * k5[i] + a76 * k6[i]);
derivs(x + c7 * h, ytemp, k7);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a81 * dydx[i] + a84 * k4[i] + a85 * k5[i] + a86 * k6[i] + a87 * k7[i]);
derivs(x + c8 * h, ytemp, k8);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a91 * dydx[i] + a94 * k4[i] + a95 * k5[i] + a96 * k6[i] + a97 * k7[i] +
a98 * k8[i]);
derivs(x + c9 * h, ytemp, k9);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a101 * dydx[i] + a104 * k4[i] + a105 * k5[i] + a106 * k6[i] +
a107 * k7[i] + a108 * k8[i] + a109 * k9[i]);
derivs(x + c10 * h, ytemp, k10);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a111 * dydx[i] + a114 * k4[i] + a115 * k5[i] + a116 * k6[i] +
a117 * k7[i] + a118 * k8[i] + a119 * k9[i] + a1110 * k10[i]);
derivs(x + c11 * h, ytemp, k2);
Doub xph = x + h;
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a121 * dydx[i] + a124 * k4[i] + a125 * k5[i] + a126 * k6[i] +
a127 * k7[i] + a128 * k8[i] + a129 * k9[i] + a1210 * k10[i] + a1211 * k2[i]);
derivs(xph, ytemp, k3);
for (i = 0; i < n; i++) {
k4[i] = b1 * dydx[i] + b6 * k6[i] + b7 * k7[i] + b8 * k8[i] + b9 * k9[i] + b10 * k10[i] +
b11 * k2[i] + b12 * k3[i];
yout[i] = y[i] + h * k4[i];
}
for (i = 0; i < n; i++) {
yerr[i] = k4[i] - bhh1 * dydx[i] - bhh2 * k9[i] - bhh3 * k3[i];
yerr2[i] = er1 * dydx[i] + er6 * k6[i] + er7 * k7[i] + er8 * k8[i] + er9 * k9[i] +
er10 * k10[i] + er11 * k2[i] + er12 * k3[i];
}
}
template<class D>
void StepperDopr853<D>::prepare_dense(const Doub h, VecDoub_I &dydxnew,
D &derivs) {
Int i;
Doub ydiff, bspl;
VecDoub ytemp(n);
for (i = 0; i < n; i++) {
rcont1[i] = y[i];
ydiff = yout[i] - y[i];
rcont2[i] = ydiff;
bspl = h * dydx[i] - ydiff;
rcont3[i] = bspl;
rcont4[i] = ydiff - h * dydxnew[i] - bspl;
rcont5[i] = d41 * dydx[i] + d46 * k6[i] + d47 * k7[i] + d48 * k8[i] +
d49 * k9[i] + d410 * k10[i] + d411 * k2[i] + d412 * k3[i];
rcont6[i] = d51 * dydx[i] + d56 * k6[i] + d57 * k7[i] + d58 * k8[i] +
d59 * k9[i] + d510 * k10[i] + d511 * k2[i] + d512 * k3[i];
rcont7[i] = d61 * dydx[i] + d66 * k6[i] + d67 * k7[i] + d68 * k8[i] +
d69 * k9[i] + d610 * k10[i] + d611 * k2[i] + d612 * k3[i];
rcont8[i] = d71 * dydx[i] + d76 * k6[i] + d77 * k7[i] + d78 * k8[i] +
d79 * k9[i] + d710 * k10[i] + d711 * k2[i] + d712 * k3[i];
}
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a141 * dydx[i] + a147 * k7[i] + a148 * k8[i] + a149 * k9[i] +
a1410 * k10[i] + a1411 * k2[i] + a1412 * k3[i] + a1413 * dydxnew[i]);
derivs(x + c14 * h, ytemp, k10);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a151 * dydx[i] + a156 * k6[i] + a157 * k7[i] + a158 * k8[i] +
a1511 * k2[i] + a1512 * k3[i] + a1513 * dydxnew[i] + a1514 * k10[i]);
derivs(x + c15 * h, ytemp, k2);
for (i = 0; i < n; i++)
ytemp[i] = y[i] + h * (a161 * dydx[i] + a166 * k6[i] + a167 * k7[i] + a168 * k8[i] +
a169 * k9[i] + a1613 * dydxnew[i] + a1614 * k10[i] + a1615 * k2[i]);
derivs(x + c16 * h, ytemp, k3);
for (i = 0; i < n; i++) {
rcont5[i] = h * (rcont5[i] + d413 * dydxnew[i] + d414 * k10[i] + d415 * k2[i] + d416 * k3[i]);
rcont6[i] = h * (rcont6[i] + d513 * dydxnew[i] + d514 * k10[i] + d515 * k2[i] + d516 * k3[i]);
rcont7[i] = h * (rcont7[i] + d613 * dydxnew[i] + d614 * k10[i] + d615 * k2[i] + d616 * k3[i]);
rcont8[i] = h * (rcont8[i] + d713 * dydxnew[i] + d714 * k10[i] + d715 * k2[i] + d716 * k3[i]);
}
}
template<class D>
Doub StepperDopr853<D>::dense_out(const Int i, const Doub x, const Doub h) {
Doub s = (x - xold) / h;
Doub s1 = 1.0 - s;
return rcont1[i] + s * (rcont2[i] + s1 * (rcont3[i] + s * (rcont4[i] + s1 * (rcont5[i] +
s * (rcont6[i] + s1 * (rcont7[i] + s *
rcont8[i]))))));
}
template<class D>
Doub StepperDopr853<D>::error(const Doub h) {
Doub err = 0.0, err2 = 0.0, sk, deno;
for (Int i = 0; i < n; i++) {
sk = atol + rtol * MAX(abs(y[i]), abs(yout[i]));
err2 += SQR(yerr[i] / sk);
err += SQR(yerr2[i] / sk);
}
deno = err + 0.01 * err2;
if (deno <= 0.0)
deno = 1.0;
return abs(h) * err * sqrt(1.0 / (n * deno));
}
template<class D>
StepperDopr853<D>::Controller::Controller() : reject(false), errold(1.0e-4) {}
template<class D>
bool StepperDopr853<D>::Controller::success(const Doub err, Doub &h) {
// static const Doub beta = 0.0, alpha = 1.0 / 8.0 - beta * 0.2, safe = 0.9, minscale = 0.333,
// maxscale = 6.0;
static const Doub beta = 0.06, alpha = 1.0 / 8.0 - beta * 0.2, safe = 0.9, minscale = 0.333,
maxscale = 6.0;
// beta = 0.04 does not really help the error (low CPU cost)
Doub scale;
if (err <= 1.0) {
if (err == 0.0)
scale = maxscale;
else {
scale = safe * pow(err, -alpha) * pow(errold, beta);
if (scale < minscale) scale = minscale;
if (scale > maxscale) scale = maxscale;
}
if (reject)
hnext = h * MIN(scale, 1.0);
else
hnext = h * scale;
errold = MAX(err, 1.0e-4);
reject = false;
return true;
} else {
scale = MAX(safe * pow(err, -alpha), minscale);
h *= scale;
reject = true;
return false;
}
}
