1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode.nonstiff;
19
20
21 import org.apache.commons.math.ode.AbstractIntegrator;
22 import org.apache.commons.math.ode.DerivativeException;
23 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
24 import org.apache.commons.math.ode.IntegratorException;
25 import org.apache.commons.math.ode.events.CombinedEventsManager;
26 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
27 import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
28 import org.apache.commons.math.ode.sampling.StepHandler;
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55 public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
56
57
58
59
60
61
62
63
64
65
66
67 protected RungeKuttaIntegrator(final String name,
68 final double[] c, final double[][] a, final double[] b,
69 final RungeKuttaStepInterpolator prototype,
70 final double step) {
71 super(name);
72 this.c = c;
73 this.a = a;
74 this.b = b;
75 this.prototype = prototype;
76 this.step = Math.abs(step);
77 }
78
79
80 public double integrate(final FirstOrderDifferentialEquations equations,
81 final double t0, final double[] y0,
82 final double t, final double[] y)
83 throws DerivativeException, IntegratorException {
84
85 sanityChecks(equations, t0, y0, t, y);
86 setEquations(equations);
87 resetEvaluations();
88 final boolean forward = (t > t0);
89
90
91 final int stages = c.length + 1;
92 if (y != y0) {
93 System.arraycopy(y0, 0, y, 0, y0.length);
94 }
95 final double[][] yDotK = new double[stages][];
96 for (int i = 0; i < stages; ++i) {
97 yDotK [i] = new double[y0.length];
98 }
99 final double[] yTmp = new double[y0.length];
100
101
102 AbstractStepInterpolator interpolator;
103 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
104 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
105 rki.reinitialize(this, yTmp, yDotK, forward);
106 interpolator = rki;
107 } else {
108 interpolator = new DummyStepInterpolator(yTmp, forward);
109 }
110 interpolator.storeTime(t0);
111
112
113 stepStart = t0;
114 stepSize = forward ? step : -step;
115 for (StepHandler handler : stepHandlers) {
116 handler.reset();
117 }
118 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
119 boolean lastStep = false;
120
121
122 while (!lastStep) {
123
124 interpolator.shift();
125
126 for (boolean loop = true; loop;) {
127
128
129 computeDerivatives(stepStart, y, yDotK[0]);
130
131
132 for (int k = 1; k < stages; ++k) {
133
134 for (int j = 0; j < y0.length; ++j) {
135 double sum = a[k-1][0] * yDotK[0][j];
136 for (int l = 1; l < k; ++l) {
137 sum += a[k-1][l] * yDotK[l][j];
138 }
139 yTmp[j] = y[j] + stepSize * sum;
140 }
141
142 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
143
144 }
145
146
147 for (int j = 0; j < y0.length; ++j) {
148 double sum = b[0] * yDotK[0][j];
149 for (int l = 1; l < stages; ++l) {
150 sum += b[l] * yDotK[l][j];
151 }
152 yTmp[j] = y[j] + stepSize * sum;
153 }
154
155
156 interpolator.storeTime(stepStart + stepSize);
157 if (manager.evaluateStep(interpolator)) {
158 final double dt = manager.getEventTime() - stepStart;
159 if (Math.abs(dt) <= Math.ulp(stepStart)) {
160
161 loop = false;
162 } else {
163
164 stepSize = dt;
165 }
166 } else {
167 loop = false;
168 }
169
170 }
171
172
173 final double nextStep = stepStart + stepSize;
174 System.arraycopy(yTmp, 0, y, 0, y0.length);
175 manager.stepAccepted(nextStep, y);
176 lastStep = manager.stop();
177
178
179 interpolator.storeTime(nextStep);
180 for (StepHandler handler : stepHandlers) {
181 handler.handleStep(interpolator, lastStep);
182 }
183 stepStart = nextStep;
184
185 if (manager.reset(stepStart, y) && ! lastStep) {
186
187
188 computeDerivatives(stepStart, y, yDotK[0]);
189 }
190
191
192 stepSize = forward ? step : -step;
193
194 }
195
196 final double stopTime = stepStart;
197 stepStart = Double.NaN;
198 stepSize = Double.NaN;
199 return stopTime;
200
201 }
202
203
204 private double[] c;
205
206
207 private double[][] a;
208
209
210 private double[] b;
211
212
213 private RungeKuttaStepInterpolator prototype;
214
215
216 private double step;
217
218 }