1 /*
2 * Licensed to the Apache Software Foundation (ASF) 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 ASF 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 * http://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.apache.commons.math4.legacy.ode;
18
19 import java.util.ArrayList;
20 import java.util.List;
21
22 import org.apache.commons.math4.legacy.core.RealFieldElement;
23 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
24 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25 import org.apache.commons.math4.legacy.core.MathArrays;
26
27
28 /**
29 * This class represents a combined set of first order differential equations,
30 * with at least a primary set of equations expandable by some sets of secondary
31 * equations.
32 * <p>
33 * One typical use case is the computation of the Jacobian matrix for some ODE.
34 * In this case, the primary set of equations corresponds to the raw ODE, and we
35 * add to this set another bunch of secondary equations which represent the Jacobian
36 * matrix of the primary set.
37 * </p>
38 * <p>
39 * We want the integrator to use <em>only</em> the primary set to estimate the
40 * errors and hence the step sizes. It should <em>not</em> use the secondary
41 * equations in this computation. The {@link FirstOrderFieldIntegrator integrator} will
42 * be able to know where the primary set ends and so where the secondary sets begin.
43 * </p>
44 *
45 * @see FirstOrderFieldDifferentialEquations
46 * @see FieldSecondaryEquations
47 *
48 * @param <T> the type of the field elements
49 * @since 3.6
50 */
51
52 public class FieldExpandableODE<T extends RealFieldElement<T>> {
53
54 /** Primary differential equation. */
55 private final FirstOrderFieldDifferentialEquations<T> primary;
56
57 /** Components of the expandable ODE. */
58 private List<FieldSecondaryEquations<T>> components;
59
60 /** Mapper for all equations. */
61 private FieldEquationsMapper<T> mapper;
62
63 /** Build an expandable set from its primary ODE set.
64 * @param primary the primary set of differential equations to be integrated.
65 */
66 public FieldExpandableODE(final FirstOrderFieldDifferentialEquations<T> primary) {
67 this.primary = primary;
68 this.components = new ArrayList<>();
69 this.mapper = new FieldEquationsMapper<>(null, primary.getDimension());
70 }
71
72 /** Get the mapper for the set of equations.
73 * @return mapper for the set of equations
74 */
75 public FieldEquationsMapper<T> getMapper() {
76 return mapper;
77 }
78
79 /** Add a set of secondary equations to be integrated along with the primary set.
80 * @param secondary secondary equations set
81 * @return index of the secondary equation in the expanded state, to be used
82 * as the parameter to {@link FieldODEState#getSecondaryState(int)} and
83 * {@link FieldODEStateAndDerivative#getSecondaryDerivative(int)} (beware index
84 * 0 corresponds to main state, additional states start at 1)
85 */
86 public int addSecondaryEquations(final FieldSecondaryEquations<T> secondary) {
87
88 components.add(secondary);
89 mapper = new FieldEquationsMapper<>(mapper, secondary.getDimension());
90
91 return components.size();
92 }
93
94 /** Initialize equations at the start of an ODE integration.
95 * @param t0 value of the independent <I>time</I> variable at integration start
96 * @param y0 array containing the value of the state vector at integration start
97 * @param finalTime target time for the integration
98 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
99 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
100 */
101 public void init(final T t0, final T[] y0, final T finalTime) {
102
103 // initialize primary equations
104 int index = 0;
105 final T[] primary0 = mapper.extractEquationData(index, y0);
106 primary.init(t0, primary0, finalTime);
107
108 // initialize secondary equations
109 while (++index < mapper.getNumberOfEquations()) {
110 final T[] secondary0 = mapper.extractEquationData(index, y0);
111 components.get(index - 1).init(t0, primary0, secondary0, finalTime);
112 }
113 }
114
115 /** Get the current time derivative of the complete state vector.
116 * @param t current value of the independent <I>time</I> variable
117 * @param y array containing the current value of the complete state vector
118 * @return time derivative of the complete state vector
119 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
120 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
121 */
122 public T[] computeDerivatives(final T t, final T[] y)
123 throws MaxCountExceededException, DimensionMismatchException {
124
125 final T[] yDot = MathArrays.buildArray(t.getField(), mapper.getTotalDimension());
126
127 // compute derivatives of the primary equations
128 int index = 0;
129 final T[] primaryState = mapper.extractEquationData(index, y);
130 final T[] primaryStateDot = primary.computeDerivatives(t, primaryState);
131 mapper.insertEquationData(index, primaryStateDot, yDot);
132
133 // Add contribution for secondary equations
134 while (++index < mapper.getNumberOfEquations()) {
135 final T[] componentState = mapper.extractEquationData(index, y);
136 final T[] componentStateDot = components.get(index - 1).computeDerivatives(t, primaryState, primaryStateDot,
137 componentState);
138 mapper.insertEquationData(index, componentStateDot, yDot);
139 }
140
141 return yDot;
142 }
143 }