37 const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
42 SIBS::redMax = 1.0e-5,
54 alpha_(kMaxX_, kMaxX_, 0.0),
55 d_p_(n_, kMaxX_, 0.0),
84 bool exitflag =
false;
88 hNext = xNew_ = -GREAT;
89 scalar eps1 = safe1*eps;
92 for (
register label
k=0;
k<kMaxX_;
k++)
94 a_[
k + 1] = a_[
k] + nSeq_[
k + 1];
97 for (
register label iq = 1; iq<kMaxX_; iq++)
99 for (
register label
k=0;
k<iq;
k++)
102 pow(eps1, (a_[
k + 1] - a_[iq + 1])
103 /((a_[iq + 1] - a_[0] + 1.0)*(2*
k + 3)));
110 for (
register label
k=0;
k<kMaxX_;
k++)
112 a_[
k + 1] = a_[
k] + nSeq_[
k + 1];
115 for (kOpt_ = 1; kOpt_<kMaxX_ - 1; kOpt_++)
117 if (a_[kOpt_ + 1] > a_[kOpt_]*alpha_[kOpt_ - 1][kOpt_])
132 if (x != xNew_ || h != hNext)
140 scalar maxErr = SMALL;
146 for (k=0; k <= kMax_; k++)
153 <<
"step size underflow"
157 SIMPR(ode, x, yTemp_, dydx, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
158 scalar xest =
sqr(h/nSeq_[k]);
160 polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
165 for (
register label i=0; i<n_; i++)
167 maxErr =
max(maxErr,
mag(yErr_[i]/yScale[i]));
171 err_[km] =
pow(maxErr/safe1, 1.0/(2*km + 3));
174 if (k != 0 && (k >= kOpt_ - 1 || first_))
182 if (k == kMax_ || k == kOpt_ + 1)
184 red = safe2/err_[km];
187 else if (k == kOpt_ && alpha_[kOpt_ - 1][kOpt_] < err_[km])
192 else if (kOpt_ == kMax_ && alpha_[km][kMax_ - 1] < err_[km])
194 red = alpha_[km][kMax_ - 1]*safe2/err_[km];
197 else if (alpha_[km][kOpt_] < err_[km])
199 red = alpha_[km][kOpt_ - 1]/err_[km];
210 red =
min(red, redMin);
211 red =
max(red, redMax);
219 scalar wrkmin = GREAT;
222 for (
register label kk=0; kk<=km; kk++)
224 scalar fact =
max(err_[kk], scaleMX);
225 scalar work = fact*a_[kk + 1];
236 if (kOpt_ >= k && kOpt_ != kMax_ && !reduct)
238 scalar fact =
max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX);
239 if (a_[kOpt_ + 1]*fact <= wrkmin)