38 KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25,
39 KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296;
42 KRR4::gamma = 1.0/2.0,
43 KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
44 KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
45 KRR4::c41 = -112.0/125.0, KRR4::c42 = -54.0/125.0, KRR4::c43 = -2.0/5.0,
46 KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0,
47 KRR4::b4 = 125.0/108.0,
48 KRR4::e1 = 17.0/54.0, KRR4::e2 = 7.0/36.0, KRR4::e3 = 0.0,
49 KRR4::e4 = 125.0/108.0,
50 KRR4::c1X = 1.0/2.0, KRR4::c2X = -3.0/2.0, KRR4::c3X = 121.0/50.0,
51 KRR4::c4X = 29.0/250.0,
52 KRR4::a2X = 1.0, KRR4::a3X = 3.0/5.0;
71 pivotIndices_(n_, 0.0)
94 ode.
jacobian(xTemp, yTemp_, dfdx_, dfdy_);
98 for (
register label jtry=0; jtry<maxtry; jtry++)
100 for (
register label i=0; i<n_; i++)
102 for (
register label j=0; j<n_; j++)
104 a_[i][j] = -dfdy_[i][j];
107 a_[i][i] += 1.0/(gamma*
h);
112 for (
register label i=0; i<n_; i++)
114 g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
119 for (
register label i=0; i<n_; i++)
121 y[i] = yTemp_[i] + a21*g1_[i];
127 for (
register label i=0; i<n_; i++)
129 g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/
h;
134 for (
register label i=0; i<n_; i++)
136 y[i] = yTemp_[i] + a31*g1_[i] + a32*g2_[i];
142 for (
register label i=0; i<n_; i++)
144 g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
149 for (
register label i=0; i<n_; i++)
151 g4_[i] = dydx_[i] + h*c4X*dfdx_[i]
152 + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
157 for (
register label i=0; i<n_; i++)
159 y[i] = yTemp_[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
160 yErr_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
168 <<
"stepsize not significant"
173 for (
register label i=0; i<n_; i++)
175 maxErr =
max(maxErr,
mag(yErr_[i]/yScale[i]));
182 hNext = (maxErr > errcon ? safety*h*
pow(maxErr, pgrow) : grow*
h);
187 hNext = safety*h*
pow(maxErr, pshrink);
188 h = (h >= 0.0 ?
max(hNext, shrink*h) :
min(hNext, shrink*h));