1 package nl.tudelft.simulation.jstats.ode;
2
3 import org.djutils.event.LocalEventProducer;
4
5 import nl.tudelft.simulation.jstats.ode.integrators.NumericalIntegrator;
6 import nl.tudelft.simulation.jstats.ode.integrators.NumericalIntegratorType;
7
8 /**
9 * The DifferentialEquation is the abstract basis for the DESS formalism.
10 * <p>
11 * Copyright (c) 2002-2024 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
12 * for project information <a href="https://simulation.tudelft.nl/" target="_blank"> https://simulation.tudelft.nl</a>. The DSOL
13 * project is distributed under a three-clause BSD-style license, which can be found at
14 * <a href="https://https://simulation.tudelft.nl/dsol/docs/latest/license.html" target="_blank">
15 * https://https://simulation.tudelft.nl/dsol/docs/latest/license.html</a>.
16 * </p>
17 * @author <a href="https://www.tudelft.nl/averbraeck" target="_blank"> Alexander Verbraeck</a>
18 * @author <a href="https://www.linkedin.com/in/peterhmjacobs">Peter Jacobs </a>
19 */
20 public abstract class DifferentialEquation extends LocalEventProducer implements DifferentialEquationInterface
21 {
22 /** */
23 private static final long serialVersionUID = 1L;
24
25 /** the numerical integrator for the differential equations. */
26 private NumericalIntegrator integrator = null;
27
28 /** the last calculated value array for lastX, initialized with the initial value array y0. */
29 @SuppressWarnings("checkstyle:visibilitymodifier")
30 protected double[] lastY = null;
31
32 /** the stepSize; can be negative or positive. */
33 @SuppressWarnings("checkstyle:visibilitymodifier")
34 protected double stepSize = Double.NaN;
35
36 /** the last x value, initialized with x0 to start integration. */
37 @SuppressWarnings("checkstyle:visibilitymodifier")
38 protected double lastX = Double.NaN;
39
40 /**
41 * constructs a new DifferentialEquation with a user-specified integrator.
42 * @param stepSize double; the stepSize to use.
43 * @param integratorType NumericalIntegratorType; the integrator to use.
44 */
45 public DifferentialEquation(final double stepSize, final NumericalIntegratorType integratorType)
46 {
47 super();
48 this.stepSize = stepSize;
49 this.integrator = integratorType.getInstance(stepSize, this);
50 }
51
52 /** {@inheritDoc} */
53 @Override
54 public void initialize(final double x0, final double[] y0)
55 {
56 this.lastX = x0;
57 this.lastY = y0;
58 }
59
60 /** {@inheritDoc} */
61 @Override
62 public double[] y(final double x)
63 {
64 // If the ODE is not initialized, the cache is empty.
65 if (Double.isNaN(this.lastX))
66 {
67 throw new RuntimeException("differential equation not initialized");
68 }
69
70 if (x == this.lastX)
71 {
72 return this.lastY;
73 }
74
75 // Are we integrating in the right direction?
76 if (Math.signum(this.stepSize) != Math.signum(x - this.lastX))
77 {
78 throw new RuntimeException("Sign of the stepsize does not integrate towards x from x0");
79 }
80 return this.integrateY(x, this.lastX, this.lastY);
81 }
82
83 /**
84 * integrates Y.
85 * @param x double; the x-value
86 * @param initialX double; the initial X value, non-final (will be updated)
87 * @param initialY double[]; the initial Y value, non-final (will be updated)
88 * @return the new Y value
89 */
90 @SuppressWarnings("checkstyle:finalparameters")
91 protected double[] integrateY(final double x, /* non-final */ double initialX, /* non-final */ double[] initialY)
92 {
93 // we request the new value from the integrator.
94 if (this.stepSize > 0)
95 {
96 while (x > initialX + this.stepSize)
97 {
98 initialY = this.integrator.next(initialX, initialY); // XXX: this process can adapt the stepSize (!)
99 initialX = initialX + this.stepSize;
100 }
101 }
102 else // negative stepsize!
103 {
104 while (x < initialX + this.stepSize)
105 {
106 initialY = this.integrator.next(initialX, initialY); // XXX: this process can adapt the stepSize (!)
107 initialX = initialX + this.stepSize;
108 }
109 }
110
111 // We are in our final step. // XXX don't understand
112 double[] nextValue = this.integrator.next(initialX, initialY);
113 double ratio = (x - initialX) / this.stepSize;
114 for (int i = 0; i < initialY.length; i++)
115 {
116 initialY[i] = initialY[i] + ratio * (nextValue[i] - initialY[i]);
117 }
118 this.lastX = x;
119 this.lastY = initialY;
120 return initialY;
121 }
122
123 /**
124 * @return Returns the integrator.
125 */
126 public NumericalIntegrator getIntegrator()
127 {
128 return this.integrator;
129 }
130
131 /**
132 * @param integrator NumericalIntegrator; The integrator to set.
133 */
134 public void setIntegrator(final NumericalIntegrator integrator)
135 {
136 this.integrator = integrator;
137 }
138
139 }