40 const int eul_max_count = 7;
41 const double eul_err_coef = 1.0/150.0;
42 const double eul_step_coef = 0.25;
43 const double eul_step_rat = 0.01;
44 const int fw_max_count = 7;
45 const double fw_err_coef = 1.0/150.0;
46 const double fw_err_rnghi = 1.5;
47 const double fw_err_rnglo = 0.75;
50 bool EulDoubleStepFlag;
51 bool FWDoubleStepFlag;
56 static int FWDoubleCount;
57 static int EulDoubleCount;
58 static double Eul_StepSize;
59 static double PrevStep;
61 Dprintf((
" Fowler-Warten integration step "));
66 FWDoubleStepFlag =
true;
67 EulDoubleStepFlag =
true;
68 FWHalveStepFlag =
false;
95 Dprintf((
"E_MIN: %g, E_MAX %g", eul_step_coef*SIMLIB_MinStep,
96 eul_step_coef*SIMLIB_StepSize));
100 (*ip)->SetState((*ip)->GetOldState()+Eul_StepSize*(*ip)->GetOldDiff());
117 eerr = Eul_StepSize*fabs((*ip)->GetDiff() - (*ip)->GetOldDiff());
120 if(eerr < eul_err_coef*terr)
123 EulDoubleStepFlag =
false;
126 if(Eul_StepSize>eul_step_coef*SIMLIB_MinStep) {
128 Eul_StepSize =
max(0.5*Eul_StepSize, eul_step_coef*SIMLIB_MinStep);
133 _Print(
"\n Integrator[%i] ",i);
155 Dprintf((
"E_F: %d, E_C: %d",EulDoubleStepFlag, EulDoubleCount));
158 if(EulDoubleStepFlag) {
165 if(EulDoubleCount >= eul_max_count) {
167 Eul_StepSize=
min(eul_step_coef*SIMLIB_StepSize, 2.0*Eul_StepSize);
170 Dprintf((
"E_S: %g", Eul_StepSize));
185 yia =
FW_First ? 0 : (((*ip)->GetOldState() -
Y1[i]) / PrevStep);
186 d1 = (*ip)->GetOldDiff() - yia;
187 d2 = ((*ip)->GetDiff() - (*ip)->GetOldDiff())/Eul_StepSize;
188 ll = (d1<=prec && d1>=-
prec) ? 0 : (d2/d1);
189 denom = SIMLIB_StepSize * ll;
190 c1 = (denom >= -
prec)
191 ? (1.0 + 0.5 * denom)
192 : ((exp(denom) - 1.0) / denom);
193 c0 = (ll>=0) ? (1.0 + denom)
196 (*ip)->SetState((*ip)->GetOldState() + SIMLIB_StepSize * (yia + c1 * d1));
197 ERR[i] = yia + c0 * d1;
214 eerr = SIMLIB_StepSize*fabs((*ip)->GetDiff() -
ERR[i]);
217 if(eerr > fw_err_rnghi*terr) {
220 FWDoubleStepFlag =
false;
221 if(SIMLIB_StepSize > SIMLIB_MinStep) {
232 _Print(
"\n Integrator[%lu] ",(
unsigned long)i);
234 else if(eerr > fw_err_rnglo*terr) {
237 FWHalveStepFlag =
true;
238 FWDoubleStepFlag =
false;
240 else if(eerr < fw_err_coef*terr) {
250 FWDoubleStepFlag = FWDoubleStepFlag && FWMayDouble;
269 Y1[i] = (*ip)->GetOldDiff();
275 if(FWHalveStepFlag) {
280 else if(FWDoubleStepFlag) {
289 if(FWDoubleCount >= fw_max_count) {
307 Dprintf((
"FW::PrepareStep()"));
320 const double FW::prec = DBL_EPSILON;
static bool IsEndStepEvent
int _Print(const char *fmt,...)
output of messages to stdout, too
void SIMLIB_Dynamic()
performs evaluation of integrators and status blocks
double SIMLIB_StepStartTime
last step time
double SIMLIB_ContractStep
requested step size
virtual bool PrepareStep(void) override
prepare object for integration step
const double & OptStep
optimal integration step
Implementation of class CalendarList interface is static - using global functions in SQS namespace...
double max(double a, double b)
double SIMLIB_RelativeError
relative error
void SIMLIB_warning(const enum _ErrEnum N)
print warning message and continue
const double & Time
model time (is NOT the block)
double SIMLIB_StepSize
actual step
virtual bool PrepareStep(void)
prepare object for integration step
static bool StateCond(void)
check on changes of state conditions
double min(double a, double b)
double SIMLIB_AbsoluteError
absolute error
bool SIMLIB_ContractStepFlag
requests shorter step
Internal header file for SIMLIB/C++.
static Iterator LastIntegrator(void)
virtual void Integrate(void) override
static Iterator FirstIntegrator(void)
Main SIMLIB/C++ interface.
IntegratorContainer::iterator Iterator
#define _SetTime(t, x)
macro for simple assignement to internal time variables
bool SIMLIB_ConditionFlag
double SIMLIB_MinStep
minimal step
double SIMLIB_DeltaTime
Time-SIMLIB_StepStartTime.
double SIMLIB_OptStep
optimal step
double SIMLIB_MaxStep
max. step