001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math3.ode; 018 019import java.util.ArrayList; 020import java.util.List; 021 022import org.apache.commons.math3.RealFieldElement; 023import org.apache.commons.math3.exception.DimensionMismatchException; 024import org.apache.commons.math3.exception.MaxCountExceededException; 025import org.apache.commons.math3.util.MathArrays; 026 027 028/** 029 * This class represents a combined set of first order differential equations, 030 * with at least a primary set of equations expandable by some sets of secondary 031 * equations. 032 * <p> 033 * One typical use case is the computation of the Jacobian matrix for some ODE. 034 * In this case, the primary set of equations corresponds to the raw ODE, and we 035 * add to this set another bunch of secondary equations which represent the Jacobian 036 * matrix of the primary set. 037 * </p> 038 * <p> 039 * We want the integrator to use <em>only</em> the primary set to estimate the 040 * errors and hence the step sizes. It should <em>not</em> use the secondary 041 * equations in this computation. The {@link FirstOrderFieldIntegrator integrator} will 042 * be able to know where the primary set ends and so where the secondary sets begin. 043 * </p> 044 * 045 * @see FirstOrderFieldDifferentialEquations 046 * @see FieldSecondaryEquations 047 * 048 * @param <T> the type of the field elements 049 * @since 3.6 050 */ 051 052public class FieldExpandableODE<T extends RealFieldElement<T>> { 053 054 /** Primary differential equation. */ 055 private final FirstOrderFieldDifferentialEquations<T> primary; 056 057 /** Components of the expandable ODE. */ 058 private List<FieldSecondaryEquations<T>> components; 059 060 /** Mapper for all equations. */ 061 private FieldEquationsMapper<T> mapper; 062 063 /** Build an expandable set from its primary ODE set. 064 * @param primary the primary set of differential equations to be integrated. 065 */ 066 public FieldExpandableODE(final FirstOrderFieldDifferentialEquations<T> primary) { 067 this.primary = primary; 068 this.components = new ArrayList<FieldSecondaryEquations<T>>(); 069 this.mapper = new FieldEquationsMapper<T>(null, primary.getDimension()); 070 } 071 072 /** Get the mapper for the set of equations. 073 * @return mapper for the set of equations 074 */ 075 public FieldEquationsMapper<T> getMapper() { 076 return mapper; 077 } 078 079 /** Add a set of secondary equations to be integrated along with the primary set. 080 * @param secondary secondary equations set 081 * @return index of the secondary equation in the expanded state, to be used 082 * as the parameter to {@link FieldODEState#getSecondaryState(int)} and 083 * {@link FieldODEStateAndDerivative#getSecondaryDerivative(int)} (beware index 084 * 0 corresponds to main state, additional states start at 1) 085 */ 086 public int addSecondaryEquations(final FieldSecondaryEquations<T> secondary) { 087 088 components.add(secondary); 089 mapper = new FieldEquationsMapper<T>(mapper, secondary.getDimension()); 090 091 return components.size(); 092 093 } 094 095 /** Initialize equations at the start of an ODE integration. 096 * @param t0 value of the independent <I>time</I> variable at integration start 097 * @param y0 array containing the value of the state vector at integration start 098 * @param finalTime target time for the integration 099 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 100 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 101 */ 102 public void init(final T t0, final T[] y0, final T finalTime) { 103 104 // initialize primary equations 105 int index = 0; 106 final T[] primary0 = mapper.extractEquationData(index, y0); 107 primary.init(t0, primary0, finalTime); 108 109 // initialize secondary equations 110 while (++index < mapper.getNumberOfEquations()) { 111 final T[] secondary0 = mapper.extractEquationData(index, y0); 112 components.get(index - 1).init(t0, primary0, secondary0, finalTime); 113 } 114 115 } 116 117 /** Get the current time derivative of the complete state vector. 118 * @param t current value of the independent <I>time</I> variable 119 * @param y array containing the current value of the complete state vector 120 * @return time derivative of the complete state vector 121 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 122 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 123 */ 124 public T[] computeDerivatives(final T t, final T[] y) 125 throws MaxCountExceededException, DimensionMismatchException { 126 127 final T[] yDot = MathArrays.buildArray(t.getField(), mapper.getTotalDimension()); 128 129 // compute derivatives of the primary equations 130 int index = 0; 131 final T[] primaryState = mapper.extractEquationData(index, y); 132 final T[] primaryStateDot = primary.computeDerivatives(t, primaryState); 133 mapper.insertEquationData(index, primaryStateDot, yDot); 134 135 // Add contribution for secondary equations 136 while (++index < mapper.getNumberOfEquations()) { 137 final T[] componentState = mapper.extractEquationData(index, y); 138 final T[] componentStateDot = components.get(index - 1).computeDerivatives(t, primaryState, primaryStateDot, 139 componentState); 140 mapper.insertEquationData(index, componentStateDot, yDot); 141 } 142 143 return yDot; 144 145 } 146 147}