1 /*
2 * Licensed to the Hipparchus project under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.hipparchus.ode;
18
19 import java.util.ArrayList;
20 import java.util.List;
21
22 import org.hipparchus.exception.MathIllegalArgumentException;
23 import org.hipparchus.exception.MathIllegalStateException;
24
25
26 /**
27 * This class represents a combined set of first order differential equations,
28 * with at least a primary set of equations expandable by some sets of secondary
29 * equations.
30 * <p>
31 * One typical use case is the computation of the Jacobian matrix for some ODE.
32 * In this case, the primary set of equations corresponds to the raw ODE, and we
33 * add to this set another bunch of secondary equations which represent the Jacobian
34 * matrix of the primary set.
35 * </p>
36 * <p>
37 * We want the integrator to use <em>only</em> the primary set to estimate the
38 * errors and hence the step sizes. It should <em>not</em> use the secondary
39 * equations in this computation. The {@link AbstractIntegrator integrator} will
40 * be able to know where the primary set ends and so where the secondary sets begin.
41 * </p>
42 *
43 * @see OrdinaryDifferentialEquation
44 * @see VariationalEquation
45 *
46 */
47
48 public class ExpandableODE {
49
50 /** Primary differential equation. */
51 private final OrdinaryDifferentialEquation primary;
52
53 /** Components of the expandable ODE. */
54 private List<SecondaryODE> components;
55
56 /** Mapper for all equations. */
57 private EquationsMapper mapper;
58
59 /** Build an expandable set from its primary ODE set.
60 * @param primary the primary set of differential equations to be integrated.
61 */
62 public ExpandableODE(final OrdinaryDifferentialEquation primary) {
63 this.primary = primary;
64 this.components = new ArrayList<>();
65 this.mapper = new EquationsMapper(null, primary.getDimension());
66 }
67
68 /** Get the primary set of differential equations to be integrated.
69 * @return primary set of differential equations to be integrated
70 */
71 public OrdinaryDifferentialEquation getPrimary() {
72 return primary;
73 }
74
75 /** Get the mapper for the set of equations.
76 * @return mapper for the set of equations
77 */
78 public EquationsMapper getMapper() {
79 return mapper;
80 }
81
82 /** Add a set of secondary equations to be integrated along with the primary set.
83 * @param secondary secondary equations set
84 * @return index of the secondary equation in the expanded state, to be used
85 * as the parameter to {@link FieldODEState#getSecondaryState(int)} and
86 * {@link FieldODEStateAndDerivative#getSecondaryDerivative(int)} (beware index
87 * 0 corresponds to primary state, secondary states start at 1)
88 */
89 public int addSecondaryEquations(final SecondaryODE secondary) {
90
91 components.add(secondary);
92 mapper = new EquationsMapper(mapper, secondary.getDimension());
93
94 return components.size();
95
96 }
97
98 /** Initialize equations at the start of an ODE integration.
99 * @param s0 state at integration start
100 * @param finalTime target time for the integration
101 * @exception MathIllegalStateException if the number of functions evaluations is exceeded
102 * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
103 */
104 public void init(final ODEState s0, final double finalTime) {
105
106 final double t0 = s0.getTime();
107
108 // initialize primary equations
109 final double[] primary0 = s0.getPrimaryState();
110 primary.init(t0, primary0, finalTime);
111
112 // initialize secondary equations
113 for (int index = 1; index < mapper.getNumberOfEquations(); ++index) {
114 final double[] secondary0 = s0.getSecondaryState(index);
115 components.get(index - 1).init(t0, primary0, secondary0, finalTime);
116 }
117
118 }
119
120 /** Get the current time derivative of the complete state vector.
121 * @param t current value of the independent <I>time</I> variable
122 * @param y array containing the current value of the complete state vector
123 * @return time derivative of the complete state vector
124 * @exception MathIllegalStateException if the number of functions evaluations is exceeded
125 * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
126 */
127 public double[] computeDerivatives(final double t, final double[] y)
128 throws MathIllegalArgumentException, MathIllegalStateException {
129
130 final double[] yDot = new double[mapper.getTotalDimension()];
131
132 // compute derivatives of the primary equations
133 final double[] primaryState = mapper.extractEquationData(0, y);
134 final double[] primaryStateDot = primary.computeDerivatives(t, primaryState);
135
136 // Add contribution for secondary equations
137 for (int index = 1; index < mapper.getNumberOfEquations(); ++index) {
138 final double[] componentState = mapper.extractEquationData(index, y);
139 final double[] componentStateDot = components.get(index - 1).computeDerivatives(t, primaryState, primaryStateDot,
140 componentState);
141 mapper.insertEquationData(index, componentStateDot, yDot);
142 }
143
144 // we retrieve the primaryStateDot array after the secondary equations have
145 // been computed in case they change the main state derivatives; this happens
146 // for example in optimal control when the secondary equations handle co-state,
147 // which changes control, and the control changes the primary state
148 mapper.insertEquationData(0, primaryStateDot, yDot);
149
150 return yDot;
151
152 }
153
154 }