|
Java example source code file (ode.xml)
The ode.xml Java example source code<?xml version="1.0"?> <!-- Licensed to the Apache Software Foundation (ASF) under one or more contributor license agreements. See the NOTICE file distributed with this work for additional information regarding copyright ownership. The ASF licenses this file to You under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. --> <?xml-stylesheet type="text/xsl" href="./xdoc.xsl"?> <document url="ode.html"> <properties> <title>The Commons Math User Guide - Ordinary Differential Equations Integration </properties> <body> <section name="15 Ordinary Differential Equations Integration"> <subsection name="15.1 Overview" href="overview"> <p> The ode package provides classes to solve Ordinary Differential Equations problems. </p> <p> This package solves Initial Value Problems of the form y'=f(t,y) with t<sub>0 and y(t<sub>0)=y0 known. The provided integrators compute an estimate of y(t) from t=t<sub>0 to t=t1. </p> <p> All integrators provide dense output. This means that besides computing the state vector at discrete times, they also provide a cheap mean to get both the state and its derivative between the time steps. They do so through classes extending the <a href="../apidocs/org/apache/commons/math3/ode/sampling/StepInterpolator.html">StepInterpolator abstract class, which are made available to the user at the end of each step. </p> <p> All integrators handle multiple discrete events detection based on switching functions. This means that the integrator can be driven by user specified discrete events (occurring when the sign of user-supplied <i>switching function changes). The steps are shortened as needed to ensure the events occur at step boundaries (even if the integrator is a fixed-step integrator). When the events are triggered, integration can be stopped (this is called a G-stop facility), the state vector can be changed, or integration can simply go on. The latter case is useful to handle discontinuities in the differential equations gracefully and get accurate dense output even close to the discontinuity. </p> <p> All integrators support setting a maximal number of evaluations of differential equations function. If this number is exceeded, an exception will be thrown during integration. This can be used to prevent infinite loops if for example error control or discrete events create a really large number of extremely small steps. By default, the maximal number of evaluation is set to <code>Integer.MAX_VALUE (i.e. 231-1 or 2147483647). It is recommended to set this maximal number to a value suited to the ODE problem, integration range, and step size or error control settings. </p> <p> All integrators support expanding the main ODE with one or more secondary ODE to manage additional state that will be integrated together with the main state. This can be used for example to integrate variational equations and compute not only the main state but also its partial derivatives with respect to either the initial state or some parameters, these derivatives being handled be secondary ODE (see below for an example). </p> <p> Two parallel APIs are available. The first is devoted to solve ode for which the integration free variable t and the state y(t) are primitive double and primitive double array respectively. Starting with version 3.6, a second API is devoted to solve ode for which the integration free variable t and the state y(t) are <code>RealFieldElement and | Name | Order | <tr>Euler | 1 | <tr>Midpoint | 2 | <tr>Classical Runge-Kutta | 4 | <tr>Gill | 4 | <tr>3/8 | 4 | <tr>Luther | 6 | </table> </p> <p> <table border="1" align="center"> <tr BGCOLOR="#CCCCFF">Adaptive Stepsize Integrators | <tr BGCOLOR="#EEEEFF">Name | Integration Order | Error Estimation Order | <tr>Higham and Hall | 5 | 4 | <tr>Dormand-Prince 5(4) | 5 | 4 | <tr>Dormand-Prince 8(5,3) | 8 | 5 and 3 | <tr>Gragg-Bulirsch-Stoer | variable (up to 18 by default) | variable | <tr>Adams-Bashforth | variable | variable | <tr>Adams-Moulton | variable | variable | </table> </p> </subsection> <subsection name="15.5 Derivatives" href="derivatives"> <p> If in addition to state y(t) the user needs to compute the sensitivity of the final state with respect to the initial state (dy/dy<sub>0) or the sensitivity of the final state with respect to some parameters of the ODE (dy/dp<sub>k), he needs to register the variational equations as a set of secondary equations appended to the main state before the integration starts. Then the integration will propagate the compound state composed of both the main state and its partial derivatives. At the end of the integration, the Jacobian matrices are extracted from the integrated secondary state. The <a href="../apidocs/org/apache/commons/math3/ode/JacobianMatrices.html">JacobianMatrices</a> class can do most of this as long as the local derivatives are provided to it. It will set up the variational equations, register them as secondary equations into the ODE, and it will set up the initial values and retrieve the intermediate and final values as Jacobian matrices. </p> <p> If for example the original state dimension is 6 and there are 3 parameters, the compound state will be a 60 elements array. The first 6 elements will be the original state, the next 36 elements will represent the 6x6 Jacobian matrix of the final state with respect to the initial state, and the remaining 18 elements will represent the 6x3 Jacobian matrix of the final state with respect to the 3 parameters. The <a href="../apidocs/org/apache/commons/math3/ode/JacobianMatrices.html">JacobianMatrices</a> class does the mapping between the 60 elements compound state and the Jacobian matrices and sets up the correcsponding secondary equations. </p> <p> As the variational equations are considered to be secondary equations here, variable step integrators ignore them for step size control: they rely only on the main state. This feature is a design choice. The rationale is to get exactly the same steps, regardless of the Jacobians being computed or not, hence ensuring reproducible results in both cases. </p> <p> What remains of user responsibility is to provide the local Jacobians df(t, y, p)/dy and df(t, y, p)/dp<sub>k corresponding the the main ODE y'=f(t, y, p). The main ODE is as usual provided by the user as a class implementing the <a href="../apidocs/org/apache/commons/math3/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations interface or a sub-interface. </p> <p> If the ODE is simple enough that the user can implement df(t, y, p)/dy directly, then instead of providing an implementation of the <a href="../apidocs/org/apache/commons/math3/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a> interface only, the user should rather provide an implementation of the <a href="../apidocs/org/apache/commons/math3/ode/MainStateJacobianProvider.html">MainStateJacobianProvider</a> interface, which extends the previous interface by adding a method to compute df(t, y, p)/dy. The user class is used as a constructor parameter of the <a href="../apidocs/org/apache/commons/math3/ode/JacobianMatrices.html">JacobianMatrices class. If the ODE is too complex or the user simply does not bother implementing df(t, y, p)/dy directly, then the ODE can still be implemented using the simple <a href="../apidocs/org/apache/commons/math3/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a> interface and given as such to another constructor of the <a href="../apidocs/org/apache/commons/math3/ode/JacobianMatrices.html">JacobianMatrices</a> class, but in this case an array hy must also be provided that will contain the step size to use form each component of the main state vector y, and the Jacobian f(t, y, p)/dy will be computed internally using finite differences. This will of course trigger more evaluations of the ODE at each step and will suffer from finite differences errors, but it is much simpler to implement from a user point of view. </p> <p> The parameters are identified by a name (a simple user defined string), which are also specified at <a href="../apidocs/org/apache/commons/math3/ode/JacobianMatrices.html">JacobianMatrices</a> class construction. If the ODE is simple enough that the user can implement df(t, y, p)/dp<sub>k directly for some of the parameters pk, then he can provide one or more classes implementing the <a href="../apidocs/org/apache/commons/math3/ode/ParameterJacobianProvider.html">ParameterJacobianProvider</a> interface by calling the JacobianMatrices.addParameterJacobianProvide method. The parameters are handled one at a time, but all the calls to ParameterJacobianProvider.computeParameterJacobian will be grouped in one sequence after the call to MainStateJacobianProvider.computeMainStateJacobian This feature can be used when all the derivatives share a lot of costly computation. In this case, the user is advised to compute all the needed derivatives at once during the call to computeMainStateJacobian, including the partial derivatives with respect to the parameters and to store the derivatives temporary. Then when the next calls to computeParameterJacobian will be triggerred, it will be sufficient to return the already computed derivatives. With this architecture, many computation can be saved. This of course implies that the classes implementing both interfaces know each other (they can even be the same class if desired, but it is not required). If the ODE is too complex or the user simply does not bother implementing df(t, y, p)/dp<sub>k directly for some k, then the JacobianMatrices.setParameterStep method should be called so finite differences are used to compute the derivatives for this parameter. It is possible to have some parameters for which derivatives are provided by a direct implementation while other parameters are computed using finite differences during the same integration. </p> <p> The following example corresponds to a simple case where all derivatives can be computed analytically. The state is a 2D point travelling along a circle. There are three parameters : the two coordinates of the center and the angular velocity. </p> <source> public class CircleODE implements MainStateJacobianProvider, ParameterJacobianProvider { public static final String CENTER_X = "cx"; public static final String CENTER_Y = "cy"; public static final String OMEGA = "omega"; private double[] c; private double omega; private double[][] savedDfDp; public CircleODE(double[] c, double omega) { this.c = c; this.omega = omega; this.savedDfDp = new double[2][3]; } public int getDimension() { return 2; } public void computeDerivatives(double t, double[] y, double[] yDot) { // the state is a 2D point, the ODE therefore corresponds to the velocity yDot[0] = omega * (c[1] - y[1]); yDot[1] = omega * (y[0] - c[0]); } public Collection<String> getParametersNames() { return Arrays.asList(CENTER_X, CENTER_Y, OMEGA); } public boolean isSupported(String name) { return CENTER_X.equals(name) || CENTER_Y.equals(name) || OMEGA.equals(name); } public void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY) { // compute the Jacobian of the main state dFdY[0][0] = 0; dFdY[0][1] = -omega; dFdY[1][0] = omega; dFdY[1][1] = 0; // precompute the derivatives with respect to the parameters, // they will be provided back when computeParameterJacobian are called later on savedDfDp[0][0] = 0; savedDfDp[0][1] = omega; savedDfDp[0][2] = c[1] - y[1]; savedDfDp[1][0] = -omega; savedDfDp[1][1] = 0; savedDfDp[1][2] = y[0] - c[0]; } public void computeParameterJacobian(double t, double[] y, double[] yDot, String paramName, double[] dFdP) { // we simply return the derivatives precomputed earlier if (CENTER_X.equals(paramName)) { dFdP[0] = savedDfDp[0][0]; dFdP[1] = savedDfDp[1][0]; } else if (CENTER_Y.equals(paramName)) { dFdP[0] = savedDfDp[0][1]; dFdP[1] = savedDfDp[1][1]; } else { dFdP[0] = savedDfDp[0][2]; dFdP[1] = savedDfDp[1][2]; } } } </source> <p> This ODE is integrated as follows: </p> <source> CircleODE circle = new CircleODE(new double[] {1.0, 1.0 }, 0.1); // here, we could select only a subset of the parameters, or use another order JacobianMatrices jm = new JacobianMatrices(circle, CircleODE.CENTER_X, CircleODE.CENTER_Y, CircleODE.OMEGA); jm.addParameterJacobianProvider(circle); ExpandableStatefulODE efode = new ExpandableStatefulODE(circle); efode.setTime(0); double[] y = { 1.0, 0.0 }; efode.setPrimaryState(y); // create the variational equations and append them to the main equations, as secondary equations jm.registerVariationalEquations(efode); // integrate the compound state, with both main and additional equations DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0e3, 1.0e-10, 1.0e-12); integrator.setMaxEvaluations(5000); integrator.integrate(efode, 20.0); // retrieve the Jacobian of the final state with respect to initial state double[][] dYdY0 = new double[2][2]; jm.getCurrentMainSetJacobian(dYdY0); // retrieve the Jacobians of the final state with resepct to the various parameters double[] dYdCx = new double[2]; double[] dYdCy = new double[2]; double[] dYdOm = new double[2]; jm.getCurrentParameterJacobian(CircleODE.CENTER_X, dYdCx); jm.getCurrentParameterJacobian(CircleODE.CENTER_Y, dYdCy); jm.getCurrentParameterJacobian(CircleODE.OMEGA, dYdOm); </source> </subsection> </section> </body> </document>
... this post is sponsored by my books ... | |
#1 New Release! |
FP Best Seller |
Copyright 1998-2024 Alvin Alexander, alvinalexander.com
All Rights Reserved.
A percentage of advertising revenue from
pages under the /java/jwarehouse
URI on this website is
paid back to open source projects.