001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.math4.legacy.ode; 019 020import java.util.ArrayList; 021import java.util.Collection; 022import java.util.Collections; 023import java.util.Comparator; 024import java.util.Iterator; 025import java.util.List; 026import java.util.SortedSet; 027import java.util.TreeSet; 028 029import org.apache.commons.math4.legacy.core.Field; 030import org.apache.commons.math4.legacy.core.RealFieldElement; 031import org.apache.commons.math4.legacy.analysis.solvers.BracketedRealFieldUnivariateSolver; 032import org.apache.commons.math4.legacy.analysis.solvers.FieldBracketingNthOrderBrentSolver; 033import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 034import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 035import org.apache.commons.math4.legacy.exception.NoBracketingException; 036import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; 037import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; 038import org.apache.commons.math4.legacy.ode.events.FieldEventHandler; 039import org.apache.commons.math4.legacy.ode.events.FieldEventState; 040import org.apache.commons.math4.legacy.ode.sampling.AbstractFieldStepInterpolator; 041import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler; 042import org.apache.commons.math4.core.jdkmath.JdkMath; 043import org.apache.commons.math4.legacy.core.IntegerSequence; 044 045/** 046 * Base class managing common boilerplate for all integrators. 047 * @param <T> the type of the field elements 048 * @since 3.6 049 */ 050public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> implements FirstOrderFieldIntegrator<T> { 051 052 /** Default relative accuracy. */ 053 private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14; 054 055 /** Default function value accuracy. */ 056 private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15; 057 058 /** Step handler. */ 059 private Collection<FieldStepHandler<T>> stepHandlers; 060 061 /** Current step start. */ 062 private FieldODEStateAndDerivative<T> stepStart; 063 064 /** Current stepsize. */ 065 private T stepSize; 066 067 /** Indicator for last step. */ 068 private boolean isLastStep; 069 070 /** Indicator that a state or derivative reset was triggered by some event. */ 071 private boolean resetOccurred; 072 073 /** Field to which the time and state vector elements belong. */ 074 private final Field<T> field; 075 076 /** Events states. */ 077 private Collection<FieldEventState<T>> eventsStates; 078 079 /** Initialization indicator of events states. */ 080 private boolean statesInitialized; 081 082 /** Name of the method. */ 083 private final String name; 084 085 /** Counter for number of evaluations. */ 086 private IntegerSequence.Incrementor evaluations; 087 088 /** Differential equations to integrate. */ 089 private transient FieldExpandableODE<T> equations; 090 091 /** Build an instance. 092 * @param field field to which the time and state vector elements belong 093 * @param name name of the method 094 */ 095 protected AbstractFieldIntegrator(final Field<T> field, final String name) { 096 this.field = field; 097 this.name = name; 098 stepHandlers = new ArrayList<>(); 099 stepStart = null; 100 stepSize = null; 101 eventsStates = new ArrayList<>(); 102 statesInitialized = false; 103 evaluations = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE); 104 } 105 106 /** Get the field to which state vector elements belong. 107 * @return field to which state vector elements belong 108 */ 109 public Field<T> getField() { 110 return field; 111 } 112 113 /** {@inheritDoc} */ 114 @Override 115 public String getName() { 116 return name; 117 } 118 119 /** {@inheritDoc} */ 120 @Override 121 public void addStepHandler(final FieldStepHandler<T> handler) { 122 stepHandlers.add(handler); 123 } 124 125 /** {@inheritDoc} */ 126 @Override 127 public Collection<FieldStepHandler<T>> getStepHandlers() { 128 return Collections.unmodifiableCollection(stepHandlers); 129 } 130 131 /** {@inheritDoc} */ 132 @Override 133 public void clearStepHandlers() { 134 stepHandlers.clear(); 135 } 136 137 /** {@inheritDoc} */ 138 @Override 139 public void addEventHandler(final FieldEventHandler<T> handler, 140 final double maxCheckInterval, 141 final double convergence, 142 final int maxIterationCount) { 143 addEventHandler(handler, maxCheckInterval, convergence, 144 maxIterationCount, 145 new FieldBracketingNthOrderBrentSolver<>(field.getZero().add(DEFAULT_RELATIVE_ACCURACY), 146 field.getZero().add(convergence), 147 field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY), 148 5)); 149 } 150 151 /** {@inheritDoc} */ 152 @Override 153 public void addEventHandler(final FieldEventHandler<T> handler, 154 final double maxCheckInterval, 155 final double convergence, 156 final int maxIterationCount, 157 final BracketedRealFieldUnivariateSolver<T> solver) { 158 eventsStates.add(new FieldEventState<>(handler, maxCheckInterval, field.getZero().add(convergence), 159 maxIterationCount, solver)); 160 } 161 162 /** {@inheritDoc} */ 163 @Override 164 public Collection<FieldEventHandler<T>> getEventHandlers() { 165 final List<FieldEventHandler<T>> list = new ArrayList<>(eventsStates.size()); 166 for (FieldEventState<T> state : eventsStates) { 167 list.add(state.getEventHandler()); 168 } 169 return Collections.unmodifiableCollection(list); 170 } 171 172 /** {@inheritDoc} */ 173 @Override 174 public void clearEventHandlers() { 175 eventsStates.clear(); 176 } 177 178 /** {@inheritDoc} */ 179 @Override 180 public FieldODEStateAndDerivative<T> getCurrentStepStart() { 181 return stepStart; 182 } 183 184 /** {@inheritDoc} */ 185 @Override 186 public T getCurrentSignedStepsize() { 187 return stepSize; 188 } 189 190 /** {@inheritDoc} */ 191 @Override 192 public void setMaxEvaluations(int maxEvaluations) { 193 evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations); 194 } 195 196 /** {@inheritDoc} */ 197 @Override 198 public int getMaxEvaluations() { 199 return evaluations.getMaximalCount(); 200 } 201 202 /** {@inheritDoc} */ 203 @Override 204 public int getEvaluations() { 205 return evaluations.getCount(); 206 } 207 208 /** Prepare the start of an integration. 209 * @param eqn equations to integrate 210 * @param t0 start value of the independent <i>time</i> variable 211 * @param y0 array containing the start value of the state vector 212 * @param t target time for the integration 213 * @return initial state with derivatives added 214 */ 215 protected FieldODEStateAndDerivative<T> initIntegration(final FieldExpandableODE<T> eqn, 216 final T t0, final T[] y0, final T t) { 217 218 this.equations = eqn; 219 evaluations = evaluations.withStart(0); 220 221 // initialize ODE 222 eqn.init(t0, y0, t); 223 224 // set up derivatives of initial state 225 final T[] y0Dot = computeDerivatives(t0, y0); 226 final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<>(t0, y0, y0Dot); 227 228 // initialize event handlers 229 for (final FieldEventState<T> state : eventsStates) { 230 state.getEventHandler().init(state0, t); 231 } 232 233 // initialize step handlers 234 for (FieldStepHandler<T> handler : stepHandlers) { 235 handler.init(state0, t); 236 } 237 238 setStateInitialized(false); 239 240 return state0; 241 } 242 243 /** Get the differential equations to integrate. 244 * @return differential equations to integrate 245 */ 246 protected FieldExpandableODE<T> getEquations() { 247 return equations; 248 } 249 250 /** Get the evaluations counter. 251 * @return evaluations counter 252 */ 253 protected IntegerSequence.Incrementor getEvaluationsCounter() { 254 return evaluations; 255 } 256 257 /** Compute the derivatives and check the number of evaluations. 258 * @param t current value of the independent <I>time</I> variable 259 * @param y array containing the current value of the state vector 260 * @return state completed with derivatives 261 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 262 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 263 * @exception NullPointerException if the ODE equations have not been set (i.e. if this method 264 * is called outside of a call to {@link #integrate(FieldExpandableODE, FieldODEState, 265 * RealFieldElement) integrate} 266 */ 267 public T[] computeDerivatives(final T t, final T[] y) 268 throws DimensionMismatchException, MaxCountExceededException, NullPointerException { 269 evaluations.increment(); 270 return equations.computeDerivatives(t, y); 271 } 272 273 /** Set the stateInitialized flag. 274 * <p>This method must be called by integrators with the value 275 * {@code false} before they start integration, so a proper lazy 276 * initialization is done automatically on the first step.</p> 277 * @param stateInitialized new value for the flag 278 */ 279 protected void setStateInitialized(final boolean stateInitialized) { 280 this.statesInitialized = stateInitialized; 281 } 282 283 /** Accept a step, triggering events and step handlers. 284 * @param interpolator step interpolator 285 * @param tEnd final integration time 286 * @return state at end of step 287 * @exception MaxCountExceededException if the interpolator throws one because 288 * the number of functions evaluations is exceeded 289 * @exception NoBracketingException if the location of an event cannot be bracketed 290 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 291 */ 292 protected FieldODEStateAndDerivative<T> acceptStep(final AbstractFieldStepInterpolator<T> interpolator, 293 final T tEnd) 294 throws MaxCountExceededException, DimensionMismatchException, NoBracketingException { 295 296 FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState(); 297 final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState(); 298 299 // initialize the events states if needed 300 if (! statesInitialized) { 301 for (FieldEventState<T> state : eventsStates) { 302 state.reinitializeBegin(interpolator); 303 } 304 statesInitialized = true; 305 } 306 307 // search for next events that may occur during the step 308 final int orderingSign = interpolator.isForward() ? +1 : -1; 309 SortedSet<FieldEventState<T>> occurringEvents = new TreeSet<>(new Comparator<FieldEventState<T>>() { 310 311 /** {@inheritDoc} */ 312 @Override 313 public int compare(FieldEventState<T> es0, FieldEventState<T> es1) { 314 return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal()); 315 } 316 }); 317 318 for (final FieldEventState<T> state : eventsStates) { 319 if (state.evaluateStep(interpolator)) { 320 // the event occurs during the current step 321 occurringEvents.add(state); 322 } 323 } 324 325 AbstractFieldStepInterpolator<T> restricted = interpolator; 326 while (!occurringEvents.isEmpty()) { 327 328 // handle the chronologically first event 329 final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator(); 330 final FieldEventState<T> currentEvent = iterator.next(); 331 iterator.remove(); 332 333 // get state at event time 334 final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime()); 335 336 // restrict the interpolator to the first part of the step, up to the event 337 restricted = restricted.restrictStep(previousState, eventState); 338 339 // advance all event states to current time 340 for (final FieldEventState<T> state : eventsStates) { 341 state.stepAccepted(eventState); 342 isLastStep = isLastStep || state.stop(); 343 } 344 345 // handle the first part of the step, up to the event 346 for (final FieldStepHandler<T> handler : stepHandlers) { 347 handler.handleStep(restricted, isLastStep); 348 } 349 350 if (isLastStep) { 351 // the event asked to stop integration 352 return eventState; 353 } 354 355 FieldODEState<T> newState = null; 356 resetOccurred = false; 357 for (final FieldEventState<T> state : eventsStates) { 358 newState = state.reset(eventState); 359 if (newState != null) { 360 // some event handler has triggered changes that 361 // invalidate the derivatives, we need to recompute them 362 final T[] y = equations.getMapper().mapState(newState); 363 final T[] yDot = computeDerivatives(newState.getTime(), y); 364 resetOccurred = true; 365 return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot); 366 } 367 } 368 369 // prepare handling of the remaining part of the step 370 previousState = eventState; 371 restricted = restricted.restrictStep(eventState, currentState); 372 373 // check if the same event occurs again in the remaining part of the step 374 if (currentEvent.evaluateStep(restricted)) { 375 // the event occurs during the current step 376 occurringEvents.add(currentEvent); 377 } 378 } 379 380 // last part of the step, after the last event 381 for (final FieldEventState<T> state : eventsStates) { 382 state.stepAccepted(currentState); 383 isLastStep = isLastStep || state.stop(); 384 } 385 isLastStep = isLastStep || currentState.getTime().subtract(tEnd).abs().getReal() <= JdkMath.ulp(tEnd.getReal()); 386 387 // handle the remaining part of the step, after all events if any 388 for (FieldStepHandler<T> handler : stepHandlers) { 389 handler.handleStep(restricted, isLastStep); 390 } 391 392 return currentState; 393 } 394 395 /** Check the integration span. 396 * @param eqn set of differential equations 397 * @param t target time for the integration 398 * @exception NumberIsTooSmallException if integration span is too small 399 * @exception DimensionMismatchException if adaptive step size integrators 400 * tolerance arrays dimensions are not compatible with equations settings 401 */ 402 protected void sanityChecks(final FieldODEState<T> eqn, final T t) 403 throws NumberIsTooSmallException, DimensionMismatchException { 404 405 final double threshold = 1000 * JdkMath.ulp(JdkMath.max(JdkMath.abs(eqn.getTime().getReal()), 406 JdkMath.abs(t.getReal()))); 407 final double dt = eqn.getTime().subtract(t).abs().getReal(); 408 if (dt <= threshold) { 409 throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL, 410 dt, threshold, false); 411 } 412 } 413 414 /** Check if a reset occurred while last step was accepted. 415 * @return true if a reset occurred while last step was accepted 416 */ 417 protected boolean resetOccurred() { 418 return resetOccurred; 419 } 420 421 /** Set the current step size. 422 * @param stepSize step size to set 423 */ 424 protected void setStepSize(final T stepSize) { 425 this.stepSize = stepSize; 426 } 427 428 /** Get the current step size. 429 * @return current step size 430 */ 431 protected T getStepSize() { 432 return stepSize; 433 } 434 /** Set current step start. 435 * @param stepStart step start 436 */ 437 protected void setStepStart(final FieldODEStateAndDerivative<T> stepStart) { 438 this.stepStart = stepStart; 439 } 440 441 /** Getcurrent step start. 442 * @return current step start 443 */ 444 protected FieldODEStateAndDerivative<T> getStepStart() { 445 return stepStart; 446 } 447 448 /** Set the last state flag. 449 * @param isLastStep if true, this step is the last one 450 */ 451 protected void setIsLastStep(final boolean isLastStep) { 452 this.isLastStep = isLastStep; 453 } 454 455 /** Check if this step is the last one. 456 * @return true if this step is the last one 457 */ 458 protected boolean isLastStep() { 459 return isLastStep; 460 } 461}