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 }