From 11fa3b6248a6e7b459442d9c82a4763115da58bf Mon Sep 17 00:00:00 2001 From: srenevey Date: Sat, 13 Jan 2024 18:41:26 -0800 Subject: [PATCH] Deployed 0f0bd66 with MkDocs version: 1.3.1 --- changelog.html | 14 ++++++++++++-- index.html | 20 +++++++------------- search/search_index.json | 2 +- sitemap.xml.gz | Bin 264 -> 264 bytes 4 files changed, 20 insertions(+), 16 deletions(-) diff --git a/changelog.html b/changelog.html index 4555250..22fce6b 100644 --- a/changelog.html +++ b/changelog.html @@ -105,6 +105,16 @@

Changelog

  • 0.3.4
    • @@ -155,7 +165,7 @@

    Changelog

  • Implement automatic stiffness detection.
  • Implement "backward" integration by allowing - to be positive or negative.
  • Fix a bug in sparse output of dopri5.
    • -
  • Add the Lorenz attractor example.
    • +
  • Add Lorenz attractor example.
    • diff --git a/index.html b/index.html index 78be2f7..877ce84 100644 --- a/index.html +++ b/index.html @@ -116,10 +116,6 @@ -
  • Acknowledgments - -
  • References @@ -136,7 +132,7 @@

    ODEs Solving in Rust

    Installation

    To start using the crate in your project, add the following dependency in your project's Cargo.toml file:

    [dependencies]
-ode_solvers = "0.3.7"
+ode_solvers = "0.4.x"

    Then, in your main file, add

    use ode_solvers::*;
@@ -145,17 +141,17 @@ 
    Type alias definition
    
 
    The numerical integration methods implemented in the crate support multi-dimensional systems. In order to define the dimension of the system, declare a type alias for the state vector. For instance
    
 
    type State = Vector3<f64>;
 
    
-
    The state representation of the system is based on the SVector<T,D> structure defined in the nalgebra crate. For convenience, ode-solvers re-exports six types to work with systems of dimension 1 to 6: Vector1<T>,..., Vector6<T>. For higher dimensions, the user should import the nalgebra crate and define a SVector<T,D>  where the second type parameter of SVector is a dimension. Note that the type T must be f64. For instance, for a 9-dimensional system, one would have:
    
+
    The state representation of the system is based on the SVector<T, D> structure defined in the nalgebra crate. For convenience, ode-solvers re-exports six types to work with systems of dimension 1 to 6: Vector1<T>,..., Vector6<T>. For higher dimensions, the user should import the nalgebra crate and define a SVector<T, D>  where the second type parameter of SVector is a dimension. For instance, for a 9-dimensional system, one would have:
    
 
    type State = SVector<f64, 9>;
 
    
 
    Alternativly, one can also use the DVector<T> structure from the nalgebra crate as the state representation. When using a DVector<T>, the number of rows in the DVector<T> defines the dimension of the system.
    
 
    type State = DVector<f64>;
 
    
 
    System definition
    
-
    The system of first order ODEs must be defined in the system method of the System<V> trait. Typically, this trait is defined for a structure containing some parameters of the model. The signature of the System<V> trait is:
    
-
    pub trait System<V> {
-    fn system(&self, x: f64, y: &V, dy: &mut V);
-    fn solout(&self, _x: f64, _y: &V, _dy: &V) -> bool {
+
    The system of first order ODEs must be defined in the system method of the System<T, V> trait. Typically, this trait is defined for a structure containing some parameters of the model. The signature of the System<T, V> trait is:
    
+
    pub trait System<T, V> {
+    fn system(&self, x: T, y: &V, dy: &mut V);
+    fn solout(&self, _x: T, _y: &V, _dy: &V) -> bool {
         false
     }
 }
@@ -280,8 +276,6 @@ 
    Adaptive Step Size
    
 
    which are combined into a single error estimator
    
 
    
 
    For a thorough discussion, see Hairer et al. [1].
    
-
    Acknowledgments
    
-
    The Dopri5 and Dop853 algorithms implemented in this crate are described in references [1] and [2]. They were originally implemented in FORTRAN by E. Hairer and G. Wanner, Université de Genève, Switzerland. This Rust implementation has been adapted from the C version written by J. Colinge, Université de Genève, Switzerland and the C++ version written by Blake Ashby, Stanford University, USA.
    
 
    References
    
 
    [1] Hairer, E., Nørsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems, Second Revised Edition, Springer, 2008
    
 
    [2] Hairer, E., Wanner, G., Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems, Second Revised Edition, Springer, 2002
    
@@ -370,5 +364,5 @@ 
 
 
diff --git a/search/search_index.json b/search/search_index.json
index c651f53..199a341 100644
--- a/search/search_index.json
+++ b/search/search_index.json
@@ -1 +1 @@
-{"config":{"indexing":"full","lang":["en"],"min_search_length":3,"prebuild_index":false,"separator":"[\\s\\-]+"},"docs":[{"location":"index.html","text":"ODEs Solving in Rust Ode-solvers is a collection of numerical methods to solve ordinary differential equations (ODEs) in Rust. The following instructions should get you up and running in no time. See the examples and the API documentation for more details. Installation To start using the crate in your project, add the following dependency in your project's Cargo.toml file: [dependencies] ode_solvers = \"0.3.7\" Then, in your main file, add use ode_solvers::*; Type alias definition The numerical integration methods implemented in the crate support multi-dimensional systems. In order to define the dimension of the system, declare a type alias for the state vector. For instance type State = Vector3; The state representation of the system is based on the SVector structure defined in the nalgebra crate. For convenience, ode-solvers re-exports six types to work with systems of dimension 1 to 6: Vector1,..., Vector6. For higher dimensions, the user should import the nalgebra crate and define a SVector where the second type parameter of SVector is a dimension. Note that the type T must be f64. For instance, for a 9-dimensional system, one would have: type State = SVector; Alternativly, one can also use the DVector structure from the nalgebra crate as the state representation. When using a DVector, the number of rows in the DVector defines the dimension of the system. type State = DVector; System definition The system of first order ODEs must be defined in the system method of the System trait. Typically, this trait is defined for a structure containing some parameters of the model. The signature of the System trait is: pub trait System { fn system(&self, x: f64, y: &V, dy: &mut V); fn solout(&self, _x: f64, _y: &V, _dy: &V) -> bool { false } } where system must contain the ODEs: the second argument is the independent variable (usually time), the third one is a vector containing the dependent variable(s), and the fourth one contains the derivative(s) of y with respect to x. The method solout is called after each successful integration step and stops the integration whenever it is evaluated as true. The implementation of that method is optional. See the examples for implementation details. Method selection The following explicit Runge-Kutta methods are implemented in the current version of the crate: Method Name Order Error estimate order Dense output order Runge-Kutta 4 Rk4 4 N/A N/A Dormand-Prince Dopri5 5 4 4 Dormand-Prince Dop853 8 (5, 3) 7 These methods are defined in the modules rk4, dopri5, and dop853. The Dormand-Prince methods feature: Adaptive step size control Automatic initial step size selection Automatic stiffness detection Sparse or dense output The first step is to bring the desired module into scope: use ode_solvers::dopri5::*; Then, a structure is created using the new or the from_param method of the corresponding struct. Refer to the API documentation for a description of the input arguments. let mut stepper = Dopri5::new(system, x0, x_end, dx, y0, rtol, atol); The system is integrated using let res = stepper.integrate(); which returns Result. Upon successful completion, res = Ok(Stats) where Stats is a structure containing some information on the integration process (such as number of function evaluations). If an error occurs, res = Err(IntegrationError). Finally, the results can be retrieved with let x_out = stepper.x_out(); let y_out = stepper.y_out(); Adaptive Step Size The adaptive step size used in Dopri5 and Dop853 is computed using a PI controller h_{n+1} = h_ns|err_n|^{-\\alpha}|err_{n-1}|^\\beta where s\u200b is a safety factor, |err_n| is the current error , and |err_{n-1}| is the previous error. The coefficients \\alpha and \\beta determine the behavior of the controller. In order to make sure that h doesn't change too abruptly, bounds on the factor of h_n are enforced: f_{min} \\leq s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{max} If s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{min}\u200b , then \u200b h_{n+1} = h_nf_{min} and if s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\geq f_{max} , h_{n+1} = h_nf_{max} . The default values of these parameters are listed in the table below. Dopri5 Dop853 \\alpha\u200b 0.17 1/8 \\beta\u200b 0.04 0 f_{min}\u200b 0.2 0.333 f_{max}\u200b 10.0 6.0 s\u200b 0.9 0.9 These default values can be overridden by using the method from_param . For Dopri5, a 5th order approximation y_1 \u200band a 4th order approximation \\hat{y}_1 are constructed using the coefficients of the method. The error estimator is then err=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} where n\u200b is the dimension of the system and sc_i is a scale factor given by sc_i=a_{tol}+r_{tol}\\max{(|y_{0i}|,|y_{1i}|)} For Dop853, an 8th order approximation y_1\u200b , a 5th order approximation \\hat{y}_1 , and a 3rd order approximation \\tilde{y}_1\u200b are constructed using the coefficients of the method. Two error estimators are computed: err_{5}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} and err_{3}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\tilde{y}_{1i}}{sc_i}} which are combined into a single error estimator err=err_5\\frac{err_5}{\\sqrt{err_5^2+0.01err_3^2}} For a thorough discussion, see Hairer et al. [1]. Acknowledgments The Dopri5 and Dop853 algorithms implemented in this crate are described in references [1] and [2]. They were originally implemented in FORTRAN by E. Hairer and G. Wanner, Universit\u00e9 de Gen\u00e8ve, Switzerland. This Rust implementation has been adapted from the C version written by J. Colinge, Universit\u00e9 de Gen\u00e8ve, Switzerland and the C++ version written by Blake Ashby, Stanford University, USA. References [1] Hairer, E., N\u00f8rsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems , Second Revised Edition, Springer, 2008 [2] Hairer, E., Wanner, G., Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems , Second Revised Edition, Springer, 2002 [3] Butcher, J.C., Numerical Methods for Ordinary Differential Equations , Third Edition, John Wiley and Sons, 2016 [4] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Explicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 17, 4 (December 1991), 533-554. [5] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Implicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 20, 4 (December 1994), 496-517. [6] S\u00f6derlind, G., Automatic Control and Adaptive Time-Stepping, Numerical Algorithms (2002) 31: 281-310.","title":"Quick reference"},{"location":"index.html#odes-solving-in-rust","text":"Ode-solvers is a collection of numerical methods to solve ordinary differential equations (ODEs) in Rust. The following instructions should get you up and running in no time. See the examples and the API documentation for more details.","title":"ODEs Solving in Rust"},{"location":"index.html#installation","text":"To start using the crate in your project, add the following dependency in your project's Cargo.toml file: [dependencies] ode_solvers = \"0.3.7\" Then, in your main file, add use ode_solvers::*;","title":"Installation"},{"location":"index.html#type-alias-definition","text":"The numerical integration methods implemented in the crate support multi-dimensional systems. In order to define the dimension of the system, declare a type alias for the state vector. For instance type State = Vector3; The state representation of the system is based on the SVector structure defined in the nalgebra crate. For convenience, ode-solvers re-exports six types to work with systems of dimension 1 to 6: Vector1,..., Vector6. For higher dimensions, the user should import the nalgebra crate and define a SVector where the second type parameter of SVector is a dimension. Note that the type T must be f64. For instance, for a 9-dimensional system, one would have: type State = SVector; Alternativly, one can also use the DVector structure from the nalgebra crate as the state representation. When using a DVector, the number of rows in the DVector defines the dimension of the system. type State = DVector;","title":"Type alias definition"},{"location":"index.html#system-definition","text":"The system of first order ODEs must be defined in the system method of the System trait. Typically, this trait is defined for a structure containing some parameters of the model. The signature of the System trait is: pub trait System { fn system(&self, x: f64, y: &V, dy: &mut V); fn solout(&self, _x: f64, _y: &V, _dy: &V) -> bool { false } } where system must contain the ODEs: the second argument is the independent variable (usually time), the third one is a vector containing the dependent variable(s), and the fourth one contains the derivative(s) of y with respect to x. The method solout is called after each successful integration step and stops the integration whenever it is evaluated as true. The implementation of that method is optional. See the examples for implementation details.","title":"System definition"},{"location":"index.html#method-selection","text":"The following explicit Runge-Kutta methods are implemented in the current version of the crate: Method Name Order Error estimate order Dense output order Runge-Kutta 4 Rk4 4 N/A N/A Dormand-Prince Dopri5 5 4 4 Dormand-Prince Dop853 8 (5, 3) 7 These methods are defined in the modules rk4, dopri5, and dop853. The Dormand-Prince methods feature: Adaptive step size control Automatic initial step size selection Automatic stiffness detection Sparse or dense output The first step is to bring the desired module into scope: use ode_solvers::dopri5::*; Then, a structure is created using the new or the from_param method of the corresponding struct. Refer to the API documentation for a description of the input arguments. let mut stepper = Dopri5::new(system, x0, x_end, dx, y0, rtol, atol); The system is integrated using let res = stepper.integrate(); which returns Result. Upon successful completion, res = Ok(Stats) where Stats is a structure containing some information on the integration process (such as number of function evaluations). If an error occurs, res = Err(IntegrationError). Finally, the results can be retrieved with let x_out = stepper.x_out(); let y_out = stepper.y_out();","title":"Method selection"},{"location":"index.html#adaptive-step-size","text":"The adaptive step size used in Dopri5 and Dop853 is computed using a PI controller h_{n+1} = h_ns|err_n|^{-\\alpha}|err_{n-1}|^\\beta where s\u200b is a safety factor, |err_n| is the current error , and |err_{n-1}| is the previous error. The coefficients \\alpha and \\beta determine the behavior of the controller. In order to make sure that h doesn't change too abruptly, bounds on the factor of h_n are enforced: f_{min} \\leq s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{max} If s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{min}\u200b , then \u200b h_{n+1} = h_nf_{min} and if s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\geq f_{max} , h_{n+1} = h_nf_{max} . The default values of these parameters are listed in the table below. Dopri5 Dop853 \\alpha\u200b 0.17 1/8 \\beta\u200b 0.04 0 f_{min}\u200b 0.2 0.333 f_{max}\u200b 10.0 6.0 s\u200b 0.9 0.9 These default values can be overridden by using the method from_param . For Dopri5, a 5th order approximation y_1 \u200band a 4th order approximation \\hat{y}_1 are constructed using the coefficients of the method. The error estimator is then err=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} where n\u200b is the dimension of the system and sc_i is a scale factor given by sc_i=a_{tol}+r_{tol}\\max{(|y_{0i}|,|y_{1i}|)} For Dop853, an 8th order approximation y_1\u200b , a 5th order approximation \\hat{y}_1 , and a 3rd order approximation \\tilde{y}_1\u200b are constructed using the coefficients of the method. Two error estimators are computed: err_{5}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} and err_{3}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\tilde{y}_{1i}}{sc_i}} which are combined into a single error estimator err=err_5\\frac{err_5}{\\sqrt{err_5^2+0.01err_3^2}} For a thorough discussion, see Hairer et al. [1].","title":"Adaptive Step Size"},{"location":"index.html#acknowledgments","text":"The Dopri5 and Dop853 algorithms implemented in this crate are described in references [1] and [2]. They were originally implemented in FORTRAN by E. Hairer and G. Wanner, Universit\u00e9 de Gen\u00e8ve, Switzerland. This Rust implementation has been adapted from the C version written by J. Colinge, Universit\u00e9 de Gen\u00e8ve, Switzerland and the C++ version written by Blake Ashby, Stanford University, USA.","title":"Acknowledgments"},{"location":"index.html#references","text":"[1] Hairer, E., N\u00f8rsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems , Second Revised Edition, Springer, 2008 [2] Hairer, E., Wanner, G., Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems , Second Revised Edition, Springer, 2002 [3] Butcher, J.C., Numerical Methods for Ordinary Differential Equations , Third Edition, John Wiley and Sons, 2016 [4] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Explicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 17, 4 (December 1991), 533-554. [5] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Implicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 20, 4 (December 1994), 496-517. [6] S\u00f6derlind, G., Automatic Control and Adaptive Time-Stepping, Numerical Algorithms (2002) 31: 281-310.","title":"References"},{"location":"changelog.html","text":"Changelog 0.3.7 Fix a bug where the final step wouldn't be computed if rejected by the controller. 0.3.6 Update the dependencies. 0.3.5 Fix a bug in stiffness detection. 0.3.4 Update to nalgebra's OVector to support DVector and SVector. Update the dependencies. 0.3.3 Fix a bug in RK4 integration method. Add unit tests. 0.3.2 Fix the iteration loop inside Runge-Kutta 4 method. 0.3.1 Implement Runge-Kutta 4 method. Update the dependencies. 0.3.0 Update to Rust 2018 edition. Implement the System trait. Implement a method to stop the integration when a condition is true. 0.2.0 Update the dependencies. Use more idiomatic Rust. 0.1.2 Change the signature of the function defining the ODE(s) to reduce the number of allocations. 0.1.1 Implement automatic stiffness detection. Implement \"backward\" integration by allowing x_{end} - x_0 to be positive or negative. Fix a bug in sparse output of dopri5. Add the Lorenz attractor example.","title":"Changelog"},{"location":"changelog.html#changelog","text":"0.3.7 Fix a bug where the final step wouldn't be computed if rejected by the controller. 0.3.6 Update the dependencies. 0.3.5 Fix a bug in stiffness detection. 0.3.4 Update to nalgebra's OVector to support DVector and SVector. Update the dependencies. 0.3.3 Fix a bug in RK4 integration method. Add unit tests. 0.3.2 Fix the iteration loop inside Runge-Kutta 4 method. 0.3.1 Implement Runge-Kutta 4 method. Update the dependencies. 0.3.0 Update to Rust 2018 edition. Implement the System trait. Implement a method to stop the integration when a condition is true. 0.2.0 Update the dependencies. Use more idiomatic Rust. 0.1.2 Change the signature of the function defining the ODE(s) to reduce the number of allocations. 0.1.1 Implement automatic stiffness detection. Implement \"backward\" integration by allowing x_{end} - x_0 to be positive or negative. Fix a bug in sparse output of dopri5. Add the Lorenz attractor example.","title":"Changelog"},{"location":"examples/kepler_orbit.html","text":"Kepler Orbit Problem Definition For this example, we consider a spacecraft on a Kepler orbit about the Earth and we want to predict its future trajectory. We assume that the elliptical orbit is described by the following classical orbital elements: Semi-major axis: a = 20,000 km Eccentricity: e = 0.7 Inclination: i = 35\u00b0 Right ascension of the ascending node: \u03a9 = 100\u00b0 Argument of perigee: \u03c9 = 65\u00b0 True anomaly: \u03bd = 30\u00b0 The period of an orbit is given by P = 2\\pi\\sqrt{\\frac{a^3}{\\mu}} where \\mu is the standard gravitational parameter. For the Earth, \u03bc = 398600.4354 km 3 /s 2 and thus P = 2.8149 \\cdot 10 4 s = 7.82 hours. The initial position of the spacecraft expressed in Cartesian coordinates is \\mathbf{r} = \\begin{bmatrix} -5007.2484 & -1444.9181 & 3628.5346 \\end{bmatrix} \\quad km and velocity \\mathbf{v} = \\begin{bmatrix} 0.7177 & -10.2241 & 0.7482 \\end{bmatrix} \\quad km/s such that the initial state vector is \\mathbf{s} = [\\mathbf{r}, \\mathbf{v}]^{T} . The equations of motion describing the evolution of the system are \\ddot{\\mathbf{r}} = -\\frac{\\mu}{r^3}\\mathbf{r} where r is the magnitude of \\mathbf{r} . Since the crate handles first order ODEs, this system must be transformed into a state space form representation. An appropriate change of variables y_1 = r_1 , y_2 = r_2 , y_3 = r_3 , y_4 = \\dot{r}_1 = v_1 , y_5 = \\dot{r}_2 = v_2 , y_6 = \\dot{r}_3 = v_3 yields \\dot{y}_1 = y_4 \\dot{y}_2 = y_5 \\dot{y}_3 = y_6 \\dot{y}_4 = -\\frac{\\mu}{r^3}y_1 \\dot{y}_5 = -\\frac{\\mu}{r^3}y_2 \\dot{y}_6 = -\\frac{\\mu}{r^3}y_3 which can be numerically integrated. Implementation The problem is solved using the Dopri5 method. We first import the crate and bring the required modules into scope: use ode_solvers::dopri5::*; use ode_solvers::*; Next, two type aliases are defined: type State = Vector6; type Time = f64; and a structure containing the parameters of the model is created: struct KeplerOrbit { mu: f64, } The System trait is then implemented for KeplerOrbit : impl ode_solvers::System for KeplerOrbit { // Equations of motion of the system fn system(&self, _t: Time, y: &State, dy: &mut State) { let r = (y[0] * y[0] + y[1] * y[1] + y[2] * y[2]).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = - self.mu * y[0] / r.powi(3); dy[4] = - self.mu * y[1] / r.powi(3); dy[5] = - self.mu * y[2] / r.powi(3); } } We now implement the main function. The system of ODEs is first created with the corresponding \\mu parameter. The orbital period is then computed and the initial state is initialized. The integration will be carried out using the 5th order Dormand-Prince method. Hence, a Dopri5 structure is created with the system structure, the initial time, the final time, the time increment between two outputs, the initial state vector, the relative tolerance, and the absolute tolerance as arguments. The integration is launched by calling the method integrate() and the result is stored in res . Finally, a check on the outcome of the integration is performed and the outputs can be extracted. fn main() { let system = KeplerOrbit {mu: 398600.435436} let a: f64 = 20000.0; let period = 2.0 * PI * (a.powi(3) / system.mu).sqrt(); let y0 = State::new(-5007.248417988539, -1444.918140151374, 3628.534606178356, 0.717716656891, -10.224093784269, 0.748229399696); let mut stepper = Dopri5::new(system, 0.0, 5.0 * period, 60.0, y0, 1.0e-10, 1.0e-10); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the output... // let path = Path::new(\"./outputs/kepler_orbit_dopri5.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } } Results The following figure shows the trajectory resulting from the integration. As expected, the orbit is closed. Some information on the integration process is provided in the stats variable: Number of function evaluations: 62835 Number of accepted steps: 10467 Number of rejected steps: 5 In order to assess the accuracy of the integration, we look at the specific orbital energy which is a constant of the motion. This specific energy is given by \\varepsilon = \\frac{v^2}{2} - \\frac{\\mu}{r} where v is the amplitude of the velocity and r is the distance from the Earth's center to the spacecraft. In theory, this value is constant as long as no external force acts on the system (which is the case here). From the intial conditions, the specific energy of this orbit is \\varepsilon_0 = -9.9650 \\; km^2/s^2 By computing \\varepsilon for each time step and plotting the difference \\varepsilon_i - \\varepsilon_0 , we obtain the following time evolution of the error: As can be seen on this figure, the specific energy is constant up to the 9th digit after the decimal point which shows that the accuracy of the integration is quite good. A higher accuracy can be achieved by selecting smaller values for rtol and atol .","title":"Kepler Orbit"},{"location":"examples/kepler_orbit.html#kepler-orbit","text":"","title":"Kepler Orbit"},{"location":"examples/kepler_orbit.html#problem-definition","text":"For this example, we consider a spacecraft on a Kepler orbit about the Earth and we want to predict its future trajectory. We assume that the elliptical orbit is described by the following classical orbital elements: Semi-major axis: a = 20,000 km Eccentricity: e = 0.7 Inclination: i = 35\u00b0 Right ascension of the ascending node: \u03a9 = 100\u00b0 Argument of perigee: \u03c9 = 65\u00b0 True anomaly: \u03bd = 30\u00b0 The period of an orbit is given by P = 2\\pi\\sqrt{\\frac{a^3}{\\mu}} where \\mu is the standard gravitational parameter. For the Earth, \u03bc = 398600.4354 km 3 /s 2 and thus P = 2.8149 \\cdot 10 4 s = 7.82 hours. The initial position of the spacecraft expressed in Cartesian coordinates is \\mathbf{r} = \\begin{bmatrix} -5007.2484 & -1444.9181 & 3628.5346 \\end{bmatrix} \\quad km and velocity \\mathbf{v} = \\begin{bmatrix} 0.7177 & -10.2241 & 0.7482 \\end{bmatrix} \\quad km/s such that the initial state vector is \\mathbf{s} = [\\mathbf{r}, \\mathbf{v}]^{T} . The equations of motion describing the evolution of the system are \\ddot{\\mathbf{r}} = -\\frac{\\mu}{r^3}\\mathbf{r} where r is the magnitude of \\mathbf{r} . Since the crate handles first order ODEs, this system must be transformed into a state space form representation. An appropriate change of variables y_1 = r_1 , y_2 = r_2 , y_3 = r_3 , y_4 = \\dot{r}_1 = v_1 , y_5 = \\dot{r}_2 = v_2 , y_6 = \\dot{r}_3 = v_3 yields \\dot{y}_1 = y_4 \\dot{y}_2 = y_5 \\dot{y}_3 = y_6 \\dot{y}_4 = -\\frac{\\mu}{r^3}y_1 \\dot{y}_5 = -\\frac{\\mu}{r^3}y_2 \\dot{y}_6 = -\\frac{\\mu}{r^3}y_3 which can be numerically integrated.","title":"Problem Definition"},{"location":"examples/kepler_orbit.html#implementation","text":"The problem is solved using the Dopri5 method. We first import the crate and bring the required modules into scope: use ode_solvers::dopri5::*; use ode_solvers::*; Next, two type aliases are defined: type State = Vector6; type Time = f64; and a structure containing the parameters of the model is created: struct KeplerOrbit { mu: f64, } The System trait is then implemented for KeplerOrbit : impl ode_solvers::System for KeplerOrbit { // Equations of motion of the system fn system(&self, _t: Time, y: &State, dy: &mut State) { let r = (y[0] * y[0] + y[1] * y[1] + y[2] * y[2]).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = - self.mu * y[0] / r.powi(3); dy[4] = - self.mu * y[1] / r.powi(3); dy[5] = - self.mu * y[2] / r.powi(3); } } We now implement the main function. The system of ODEs is first created with the corresponding \\mu parameter. The orbital period is then computed and the initial state is initialized. The integration will be carried out using the 5th order Dormand-Prince method. Hence, a Dopri5 structure is created with the system structure, the initial time, the final time, the time increment between two outputs, the initial state vector, the relative tolerance, and the absolute tolerance as arguments. The integration is launched by calling the method integrate() and the result is stored in res . Finally, a check on the outcome of the integration is performed and the outputs can be extracted. fn main() { let system = KeplerOrbit {mu: 398600.435436} let a: f64 = 20000.0; let period = 2.0 * PI * (a.powi(3) / system.mu).sqrt(); let y0 = State::new(-5007.248417988539, -1444.918140151374, 3628.534606178356, 0.717716656891, -10.224093784269, 0.748229399696); let mut stepper = Dopri5::new(system, 0.0, 5.0 * period, 60.0, y0, 1.0e-10, 1.0e-10); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the output... // let path = Path::new(\"./outputs/kepler_orbit_dopri5.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } }","title":"Implementation"},{"location":"examples/kepler_orbit.html#results","text":"The following figure shows the trajectory resulting from the integration. As expected, the orbit is closed. Some information on the integration process is provided in the stats variable: Number of function evaluations: 62835 Number of accepted steps: 10467 Number of rejected steps: 5 In order to assess the accuracy of the integration, we look at the specific orbital energy which is a constant of the motion. This specific energy is given by \\varepsilon = \\frac{v^2}{2} - \\frac{\\mu}{r} where v is the amplitude of the velocity and r is the distance from the Earth's center to the spacecraft. In theory, this value is constant as long as no external force acts on the system (which is the case here). From the intial conditions, the specific energy of this orbit is \\varepsilon_0 = -9.9650 \\; km^2/s^2 By computing \\varepsilon for each time step and plotting the difference \\varepsilon_i - \\varepsilon_0 , we obtain the following time evolution of the error: As can be seen on this figure, the specific energy is constant up to the 9th digit after the decimal point which shows that the accuracy of the integration is quite good. A higher accuracy can be achieved by selecting smaller values for rtol and atol .","title":"Results"},{"location":"examples/three_body_problem.html","text":"Three-Body Problem Problem Definition In this problem, we consider the circular restricted three body problem in which a spacecraft evolves under the influence of the Earth and the Moon. The second order nondimensional equations of motion of this problem are \\ddot{x}-2\\dot{y}-x = -\\frac{(1-\\mu)(x+\\mu)}{d^3}-\\frac{\\mu(x-1+\\mu)}{r^3} \\ddot{y}+2\\dot{x}-y = -\\frac{(1-\\mu)y}{d^3}-\\frac{\\mu y}{r^3} \\ddot{z} = -\\frac{(1-\\mu)z}{d^3}-\\frac{\\mu z}{r^3} where d=\\left[(x+\\mu)^2+y^2+z^2\\right]^{\\frac{1}{2}} r=\\left[(x-1+\\mu)^2+y^2+z^2\\right] \\mu is a parameter of the system. For the Earth-Moon system, \\mu=0.0123 . Rewriting this system in state space form yields \\dot{x}_1=x_4 \\dot{x}_2=x_5 \\dot{x}_3=x_6 \\dot{x}_4=x_1+2x_5-\\frac{(1-\\mu)(x_1+\\mu)}{d^3}-\\frac{\\mu(x_1-1+\\mu)}{r^3} \\dot{x}_5=-2x_4+x_2-\\frac{(1-\\mu)x_2}{d^3}-\\frac{\\mu x_2}{r^3} \\dot{x}_5=-\\frac{(1-\\mu)x_3}{d^3}-\\frac{\\mu x_3}{r^3} Implementation The problem is solved using Dop853. The crate is first imported and the modules are brought into scope: use ode_solvers::dop853::*; use ode_solvers::*; Type aliases are then defined for the state vector and the time: type State = Vector6; type Time = f64; Next, we define the problem specific structure which will contain a single parameter: struct ThreeBodyProblem { mu: f64, } for which we implement the System trait: impl ode_solvers::System for ThreeBodyProblem { fn system(&self, _t: Time, y: &State, dy: &mut State) { let d = ((y[0] + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); let r = ((y[0] - 1.0 + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = y[0] + 2.0 * y[4] - (1.0 - self.mu) * (y[0] + self.mu) / d.powi(3) - self.mu * (y[0] - 1.0 + self.mu) / r.powi(3); dy[4] = -2.0 * y[3] + y[1] - (1.0 - self.mu) * y[1] / d.powi(3) - self.mu * y[1] / r.powi(3); dy[5] = -(1.0 - self.mu) * y[2] / d.powi(3) - self.mu * y[2] / r.powi(3); } } The system is created together with the intial state vector of the spacecraft in the main function. A Dop853 structure is created and the integrate() function is called. The result is handled with a match statement. fn main() { let system = ThreeBodyProblem {mu: 0.012300118882173}; let y0 = State::new(-0.271, -0.42, 0.0, 0.3, -1.0, 0.0); let mut stepper = Dop853::new(system, 0.0, 150.0, 0.001, y0, 1.0e-14, 1.0e-14); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the results... // let path = Path::new(\"./outputs/three_body_dop853.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } } Results The trajectory of the spacecraft for a duration of 150 (dimensionless time) is shown on the following figure: On that figure, the Earth is located at (0, 0) and the Moon is at (0.9877, 0). We see that the spacecraft is first on an orbit about the Earth before moving to the vicinity of the Moon. The accuracy of the integration can be assessed by looking at the Jacobi constant. The Jacobi constant is given by C=\\frac{2(1-\\mu)}{d} + \\frac{2\\mu}{r} + x^2 +y^2 - v^2 and is an integral of the motion (i.e. is constant as long as no external force is applied). Based on the initial conditions, the Jacobi constant yields C_0=3.1830 By computing the Jacobi constant at each time step and plotting the difference C-C_0 , we get the figure below: As can be seen on that figure, the Jacobi constant is constant up to the 11th digit after the decimal point. This accuracy is quite reasonable for that kind of problem but if a higher accuracy is required, the relative tolerance and/or the absolute tolerance could be reduced.","title":"Three-body Problem"},{"location":"examples/three_body_problem.html#three-body-problem","text":"","title":"Three-Body Problem"},{"location":"examples/three_body_problem.html#problem-definition","text":"In this problem, we consider the circular restricted three body problem in which a spacecraft evolves under the influence of the Earth and the Moon. The second order nondimensional equations of motion of this problem are \\ddot{x}-2\\dot{y}-x = -\\frac{(1-\\mu)(x+\\mu)}{d^3}-\\frac{\\mu(x-1+\\mu)}{r^3} \\ddot{y}+2\\dot{x}-y = -\\frac{(1-\\mu)y}{d^3}-\\frac{\\mu y}{r^3} \\ddot{z} = -\\frac{(1-\\mu)z}{d^3}-\\frac{\\mu z}{r^3} where d=\\left[(x+\\mu)^2+y^2+z^2\\right]^{\\frac{1}{2}} r=\\left[(x-1+\\mu)^2+y^2+z^2\\right] \\mu is a parameter of the system. For the Earth-Moon system, \\mu=0.0123 . Rewriting this system in state space form yields \\dot{x}_1=x_4 \\dot{x}_2=x_5 \\dot{x}_3=x_6 \\dot{x}_4=x_1+2x_5-\\frac{(1-\\mu)(x_1+\\mu)}{d^3}-\\frac{\\mu(x_1-1+\\mu)}{r^3} \\dot{x}_5=-2x_4+x_2-\\frac{(1-\\mu)x_2}{d^3}-\\frac{\\mu x_2}{r^3} \\dot{x}_5=-\\frac{(1-\\mu)x_3}{d^3}-\\frac{\\mu x_3}{r^3}","title":"Problem Definition"},{"location":"examples/three_body_problem.html#implementation","text":"The problem is solved using Dop853. The crate is first imported and the modules are brought into scope: use ode_solvers::dop853::*; use ode_solvers::*; Type aliases are then defined for the state vector and the time: type State = Vector6; type Time = f64; Next, we define the problem specific structure which will contain a single parameter: struct ThreeBodyProblem { mu: f64, } for which we implement the System trait: impl ode_solvers::System for ThreeBodyProblem { fn system(&self, _t: Time, y: &State, dy: &mut State) { let d = ((y[0] + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); let r = ((y[0] - 1.0 + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = y[0] + 2.0 * y[4] - (1.0 - self.mu) * (y[0] + self.mu) / d.powi(3) - self.mu * (y[0] - 1.0 + self.mu) / r.powi(3); dy[4] = -2.0 * y[3] + y[1] - (1.0 - self.mu) * y[1] / d.powi(3) - self.mu * y[1] / r.powi(3); dy[5] = -(1.0 - self.mu) * y[2] / d.powi(3) - self.mu * y[2] / r.powi(3); } } The system is created together with the intial state vector of the spacecraft in the main function. A Dop853 structure is created and the integrate() function is called. The result is handled with a match statement. fn main() { let system = ThreeBodyProblem {mu: 0.012300118882173}; let y0 = State::new(-0.271, -0.42, 0.0, 0.3, -1.0, 0.0); let mut stepper = Dop853::new(system, 0.0, 150.0, 0.001, y0, 1.0e-14, 1.0e-14); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the results... // let path = Path::new(\"./outputs/three_body_dop853.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } }","title":"Implementation"},{"location":"examples/three_body_problem.html#results","text":"The trajectory of the spacecraft for a duration of 150 (dimensionless time) is shown on the following figure: On that figure, the Earth is located at (0, 0) and the Moon is at (0.9877, 0). We see that the spacecraft is first on an orbit about the Earth before moving to the vicinity of the Moon. The accuracy of the integration can be assessed by looking at the Jacobi constant. The Jacobi constant is given by C=\\frac{2(1-\\mu)}{d} + \\frac{2\\mu}{r} + x^2 +y^2 - v^2 and is an integral of the motion (i.e. is constant as long as no external force is applied). Based on the initial conditions, the Jacobi constant yields C_0=3.1830 By computing the Jacobi constant at each time step and plotting the difference C-C_0 , we get the figure below: As can be seen on that figure, the Jacobi constant is constant up to the 11th digit after the decimal point. This accuracy is quite reasonable for that kind of problem but if a higher accuracy is required, the relative tolerance and/or the absolute tolerance could be reduced.","title":"Results"}]}
\ No newline at end of file
+{"config":{"indexing":"full","lang":["en"],"min_search_length":3,"prebuild_index":false,"separator":"[\\s\\-]+"},"docs":[{"location":"index.html","text":"ODEs Solving in Rust Ode-solvers is a collection of numerical methods to solve ordinary differential equations (ODEs) in Rust. The following instructions should get you up and running in no time. See the examples and the API documentation for more details. Installation To start using the crate in your project, add the following dependency in your project's Cargo.toml file: [dependencies] ode_solvers = \"0.4.x\" Then, in your main file, add use ode_solvers::*; Type alias definition The numerical integration methods implemented in the crate support multi-dimensional systems. In order to define the dimension of the system, declare a type alias for the state vector. For instance type State = Vector3; The state representation of the system is based on the SVector structure defined in the nalgebra crate. For convenience, ode-solvers re-exports six types to work with systems of dimension 1 to 6: Vector1,..., Vector6. For higher dimensions, the user should import the nalgebra crate and define a SVector where the second type parameter of SVector is a dimension. For instance, for a 9-dimensional system, one would have: type State = SVector; Alternativly, one can also use the DVector structure from the nalgebra crate as the state representation. When using a DVector, the number of rows in the DVector defines the dimension of the system. type State = DVector; System definition The system of first order ODEs must be defined in the system method of the System trait. Typically, this trait is defined for a structure containing some parameters of the model. The signature of the System trait is: pub trait System { fn system(&self, x: T, y: &V, dy: &mut V); fn solout(&self, _x: T, _y: &V, _dy: &V) -> bool { false } } where system must contain the ODEs: the second argument is the independent variable (usually time), the third one is a vector containing the dependent variable(s), and the fourth one contains the derivative(s) of y with respect to x. The method solout is called after each successful integration step and stops the integration whenever it is evaluated as true. The implementation of that method is optional. See the examples for implementation details. Method selection The following explicit Runge-Kutta methods are implemented in the current version of the crate: Method Name Order Error estimate order Dense output order Runge-Kutta 4 Rk4 4 N/A N/A Dormand-Prince Dopri5 5 4 4 Dormand-Prince Dop853 8 (5, 3) 7 These methods are defined in the modules rk4, dopri5, and dop853. The Dormand-Prince methods feature: Adaptive step size control Automatic initial step size selection Automatic stiffness detection Sparse or dense output The first step is to bring the desired module into scope: use ode_solvers::dopri5::*; Then, a structure is created using the new or the from_param method of the corresponding struct. Refer to the API documentation for a description of the input arguments. let mut stepper = Dopri5::new(system, x0, x_end, dx, y0, rtol, atol); The system is integrated using let res = stepper.integrate(); which returns Result. Upon successful completion, res = Ok(Stats) where Stats is a structure containing some information on the integration process (such as number of function evaluations). If an error occurs, res = Err(IntegrationError). Finally, the results can be retrieved with let x_out = stepper.x_out(); let y_out = stepper.y_out(); Adaptive Step Size The adaptive step size used in Dopri5 and Dop853 is computed using a PI controller h_{n+1} = h_ns|err_n|^{-\\alpha}|err_{n-1}|^\\beta where s\u200b is a safety factor, |err_n| is the current error , and |err_{n-1}| is the previous error. The coefficients \\alpha and \\beta determine the behavior of the controller. In order to make sure that h doesn't change too abruptly, bounds on the factor of h_n are enforced: f_{min} \\leq s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{max} If s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{min}\u200b , then \u200b h_{n+1} = h_nf_{min} and if s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\geq f_{max} , h_{n+1} = h_nf_{max} . The default values of these parameters are listed in the table below. Dopri5 Dop853 \\alpha\u200b 0.17 1/8 \\beta\u200b 0.04 0 f_{min}\u200b 0.2 0.333 f_{max}\u200b 10.0 6.0 s\u200b 0.9 0.9 These default values can be overridden by using the method from_param . For Dopri5, a 5th order approximation y_1 \u200band a 4th order approximation \\hat{y}_1 are constructed using the coefficients of the method. The error estimator is then err=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} where n\u200b is the dimension of the system and sc_i is a scale factor given by sc_i=a_{tol}+r_{tol}\\max{(|y_{0i}|,|y_{1i}|)} For Dop853, an 8th order approximation y_1\u200b , a 5th order approximation \\hat{y}_1 , and a 3rd order approximation \\tilde{y}_1\u200b are constructed using the coefficients of the method. Two error estimators are computed: err_{5}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} and err_{3}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\tilde{y}_{1i}}{sc_i}} which are combined into a single error estimator err=err_5\\frac{err_5}{\\sqrt{err_5^2+0.01err_3^2}} For a thorough discussion, see Hairer et al. [1]. References [1] Hairer, E., N\u00f8rsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems , Second Revised Edition, Springer, 2008 [2] Hairer, E., Wanner, G., Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems , Second Revised Edition, Springer, 2002 [3] Butcher, J.C., Numerical Methods for Ordinary Differential Equations , Third Edition, John Wiley and Sons, 2016 [4] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Explicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 17, 4 (December 1991), 533-554. [5] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Implicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 20, 4 (December 1994), 496-517. [6] S\u00f6derlind, G., Automatic Control and Adaptive Time-Stepping, Numerical Algorithms (2002) 31: 281-310.","title":"Quick reference"},{"location":"index.html#odes-solving-in-rust","text":"Ode-solvers is a collection of numerical methods to solve ordinary differential equations (ODEs) in Rust. The following instructions should get you up and running in no time. See the examples and the API documentation for more details.","title":"ODEs Solving in Rust"},{"location":"index.html#installation","text":"To start using the crate in your project, add the following dependency in your project's Cargo.toml file: [dependencies] ode_solvers = \"0.4.x\" Then, in your main file, add use ode_solvers::*;","title":"Installation"},{"location":"index.html#type-alias-definition","text":"The numerical integration methods implemented in the crate support multi-dimensional systems. In order to define the dimension of the system, declare a type alias for the state vector. For instance type State = Vector3; The state representation of the system is based on the SVector structure defined in the nalgebra crate. For convenience, ode-solvers re-exports six types to work with systems of dimension 1 to 6: Vector1,..., Vector6. For higher dimensions, the user should import the nalgebra crate and define a SVector where the second type parameter of SVector is a dimension. For instance, for a 9-dimensional system, one would have: type State = SVector; Alternativly, one can also use the DVector structure from the nalgebra crate as the state representation. When using a DVector, the number of rows in the DVector defines the dimension of the system. type State = DVector;","title":"Type alias definition"},{"location":"index.html#system-definition","text":"The system of first order ODEs must be defined in the system method of the System trait. Typically, this trait is defined for a structure containing some parameters of the model. The signature of the System trait is: pub trait System { fn system(&self, x: T, y: &V, dy: &mut V); fn solout(&self, _x: T, _y: &V, _dy: &V) -> bool { false } } where system must contain the ODEs: the second argument is the independent variable (usually time), the third one is a vector containing the dependent variable(s), and the fourth one contains the derivative(s) of y with respect to x. The method solout is called after each successful integration step and stops the integration whenever it is evaluated as true. The implementation of that method is optional. See the examples for implementation details.","title":"System definition"},{"location":"index.html#method-selection","text":"The following explicit Runge-Kutta methods are implemented in the current version of the crate: Method Name Order Error estimate order Dense output order Runge-Kutta 4 Rk4 4 N/A N/A Dormand-Prince Dopri5 5 4 4 Dormand-Prince Dop853 8 (5, 3) 7 These methods are defined in the modules rk4, dopri5, and dop853. The Dormand-Prince methods feature: Adaptive step size control Automatic initial step size selection Automatic stiffness detection Sparse or dense output The first step is to bring the desired module into scope: use ode_solvers::dopri5::*; Then, a structure is created using the new or the from_param method of the corresponding struct. Refer to the API documentation for a description of the input arguments. let mut stepper = Dopri5::new(system, x0, x_end, dx, y0, rtol, atol); The system is integrated using let res = stepper.integrate(); which returns Result. Upon successful completion, res = Ok(Stats) where Stats is a structure containing some information on the integration process (such as number of function evaluations). If an error occurs, res = Err(IntegrationError). Finally, the results can be retrieved with let x_out = stepper.x_out(); let y_out = stepper.y_out();","title":"Method selection"},{"location":"index.html#adaptive-step-size","text":"The adaptive step size used in Dopri5 and Dop853 is computed using a PI controller h_{n+1} = h_ns|err_n|^{-\\alpha}|err_{n-1}|^\\beta where s\u200b is a safety factor, |err_n| is the current error , and |err_{n-1}| is the previous error. The coefficients \\alpha and \\beta determine the behavior of the controller. In order to make sure that h doesn't change too abruptly, bounds on the factor of h_n are enforced: f_{min} \\leq s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{max} If s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\leq f_{min}\u200b , then \u200b h_{n+1} = h_nf_{min} and if s|err_n|^{-\\alpha}|err_{n-1}|^\\beta \\geq f_{max} , h_{n+1} = h_nf_{max} . The default values of these parameters are listed in the table below. Dopri5 Dop853 \\alpha\u200b 0.17 1/8 \\beta\u200b 0.04 0 f_{min}\u200b 0.2 0.333 f_{max}\u200b 10.0 6.0 s\u200b 0.9 0.9 These default values can be overridden by using the method from_param . For Dopri5, a 5th order approximation y_1 \u200band a 4th order approximation \\hat{y}_1 are constructed using the coefficients of the method. The error estimator is then err=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} where n\u200b is the dimension of the system and sc_i is a scale factor given by sc_i=a_{tol}+r_{tol}\\max{(|y_{0i}|,|y_{1i}|)} For Dop853, an 8th order approximation y_1\u200b , a 5th order approximation \\hat{y}_1 , and a 3rd order approximation \\tilde{y}_1\u200b are constructed using the coefficients of the method. Two error estimators are computed: err_{5}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\hat{y}_{1i}}{sc_i}} and err_{3}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\frac{y_{1i}-\\tilde{y}_{1i}}{sc_i}} which are combined into a single error estimator err=err_5\\frac{err_5}{\\sqrt{err_5^2+0.01err_3^2}} For a thorough discussion, see Hairer et al. [1].","title":"Adaptive Step Size"},{"location":"index.html#references","text":"[1] Hairer, E., N\u00f8rsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems , Second Revised Edition, Springer, 2008 [2] Hairer, E., Wanner, G., Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems , Second Revised Edition, Springer, 2002 [3] Butcher, J.C., Numerical Methods for Ordinary Differential Equations , Third Edition, John Wiley and Sons, 2016 [4] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Explicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 17, 4 (December 1991), 533-554. [5] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Implicit Runge-Kutta Methods. ACM Transactions on Mathematical Software , 20, 4 (December 1994), 496-517. [6] S\u00f6derlind, G., Automatic Control and Adaptive Time-Stepping, Numerical Algorithms (2002) 31: 281-310.","title":"References"},{"location":"changelog.html","text":"Changelog 0.4.0 Add generic type for the independent variable. Update internal representation of the solver result. Add bouncing ball example. 0.3.8 Update the Butcher tableaus to declare the coefficients as constant instead of heap allocating. 0.3.7 Fix a bug where the final step wouldn't be computed if rejected by the controller. 0.3.6 Update the dependencies. 0.3.5 Fix a bug in stiffness detection. 0.3.4 Update the methods to use nalgebra's OVector to support DVector and SVector. Update the dependencies. 0.3.3 Fix a bug in RK4 integration method. Add unit tests. 0.3.2 Fix the iteration loop inside Runge-Kutta 4 method. 0.3.1 Implement Runge-Kutta 4 method. Update the dependencies. 0.3.0 Update to Rust 2018 edition. Implement the System trait. Implement a method to stop the integration when a condition is true. 0.2.0 Update the dependencies. Use more idiomatic Rust. 0.1.2 Change the signature of the function defining the ODE(s) to reduce the number of allocations. 0.1.1 Implement automatic stiffness detection. Implement \"backward\" integration by allowing x_{end} - x_0 to be positive or negative. Fix a bug in sparse output of dopri5. Add Lorenz attractor example.","title":"Changelog"},{"location":"changelog.html#changelog","text":"0.4.0 Add generic type for the independent variable. Update internal representation of the solver result. Add bouncing ball example. 0.3.8 Update the Butcher tableaus to declare the coefficients as constant instead of heap allocating. 0.3.7 Fix a bug where the final step wouldn't be computed if rejected by the controller. 0.3.6 Update the dependencies. 0.3.5 Fix a bug in stiffness detection. 0.3.4 Update the methods to use nalgebra's OVector to support DVector and SVector. Update the dependencies. 0.3.3 Fix a bug in RK4 integration method. Add unit tests. 0.3.2 Fix the iteration loop inside Runge-Kutta 4 method. 0.3.1 Implement Runge-Kutta 4 method. Update the dependencies. 0.3.0 Update to Rust 2018 edition. Implement the System trait. Implement a method to stop the integration when a condition is true. 0.2.0 Update the dependencies. Use more idiomatic Rust. 0.1.2 Change the signature of the function defining the ODE(s) to reduce the number of allocations. 0.1.1 Implement automatic stiffness detection. Implement \"backward\" integration by allowing x_{end} - x_0 to be positive or negative. Fix a bug in sparse output of dopri5. Add Lorenz attractor example.","title":"Changelog"},{"location":"examples/kepler_orbit.html","text":"Kepler Orbit Problem Definition For this example, we consider a spacecraft on a Kepler orbit about the Earth and we want to predict its future trajectory. We assume that the elliptical orbit is described by the following classical orbital elements: Semi-major axis: a = 20,000 km Eccentricity: e = 0.7 Inclination: i = 35\u00b0 Right ascension of the ascending node: \u03a9 = 100\u00b0 Argument of perigee: \u03c9 = 65\u00b0 True anomaly: \u03bd = 30\u00b0 The period of an orbit is given by P = 2\\pi\\sqrt{\\frac{a^3}{\\mu}} where \\mu is the standard gravitational parameter. For the Earth, \u03bc = 398600.4354 km 3 /s 2 and thus P = 2.8149 \\cdot 10 4 s = 7.82 hours. The initial position of the spacecraft expressed in Cartesian coordinates is \\mathbf{r} = \\begin{bmatrix} -5007.2484 & -1444.9181 & 3628.5346 \\end{bmatrix} \\quad km and velocity \\mathbf{v} = \\begin{bmatrix} 0.7177 & -10.2241 & 0.7482 \\end{bmatrix} \\quad km/s such that the initial state vector is \\mathbf{s} = [\\mathbf{r}, \\mathbf{v}]^{T} . The equations of motion describing the evolution of the system are \\ddot{\\mathbf{r}} = -\\frac{\\mu}{r^3}\\mathbf{r} where r is the magnitude of \\mathbf{r} . Since the crate handles first order ODEs, this system must be transformed into a state space form representation. An appropriate change of variables y_1 = r_1 , y_2 = r_2 , y_3 = r_3 , y_4 = \\dot{r}_1 = v_1 , y_5 = \\dot{r}_2 = v_2 , y_6 = \\dot{r}_3 = v_3 yields \\dot{y}_1 = y_4 \\dot{y}_2 = y_5 \\dot{y}_3 = y_6 \\dot{y}_4 = -\\frac{\\mu}{r^3}y_1 \\dot{y}_5 = -\\frac{\\mu}{r^3}y_2 \\dot{y}_6 = -\\frac{\\mu}{r^3}y_3 which can be numerically integrated. Implementation The problem is solved using the Dopri5 method. We first import the crate and bring the required modules into scope: use ode_solvers::dopri5::*; use ode_solvers::*; Next, two type aliases are defined: type State = Vector6; type Time = f64; and a structure containing the parameters of the model is created: struct KeplerOrbit { mu: f64, } The System trait is then implemented for KeplerOrbit : impl ode_solvers::System for KeplerOrbit { // Equations of motion of the system fn system(&self, _t: Time, y: &State, dy: &mut State) { let r = (y[0] * y[0] + y[1] * y[1] + y[2] * y[2]).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = - self.mu * y[0] / r.powi(3); dy[4] = - self.mu * y[1] / r.powi(3); dy[5] = - self.mu * y[2] / r.powi(3); } } We now implement the main function. The system of ODEs is first created with the corresponding \\mu parameter. The orbital period is then computed and the initial state is initialized. The integration will be carried out using the 5th order Dormand-Prince method. Hence, a Dopri5 structure is created with the system structure, the initial time, the final time, the time increment between two outputs, the initial state vector, the relative tolerance, and the absolute tolerance as arguments. The integration is launched by calling the method integrate() and the result is stored in res . Finally, a check on the outcome of the integration is performed and the outputs can be extracted. fn main() { let system = KeplerOrbit {mu: 398600.435436} let a: f64 = 20000.0; let period = 2.0 * PI * (a.powi(3) / system.mu).sqrt(); let y0 = State::new(-5007.248417988539, -1444.918140151374, 3628.534606178356, 0.717716656891, -10.224093784269, 0.748229399696); let mut stepper = Dopri5::new(system, 0.0, 5.0 * period, 60.0, y0, 1.0e-10, 1.0e-10); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the output... // let path = Path::new(\"./outputs/kepler_orbit_dopri5.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } } Results The following figure shows the trajectory resulting from the integration. As expected, the orbit is closed. Some information on the integration process is provided in the stats variable: Number of function evaluations: 62835 Number of accepted steps: 10467 Number of rejected steps: 5 In order to assess the accuracy of the integration, we look at the specific orbital energy which is a constant of the motion. This specific energy is given by \\varepsilon = \\frac{v^2}{2} - \\frac{\\mu}{r} where v is the amplitude of the velocity and r is the distance from the Earth's center to the spacecraft. In theory, this value is constant as long as no external force acts on the system (which is the case here). From the intial conditions, the specific energy of this orbit is \\varepsilon_0 = -9.9650 \\; km^2/s^2 By computing \\varepsilon for each time step and plotting the difference \\varepsilon_i - \\varepsilon_0 , we obtain the following time evolution of the error: As can be seen on this figure, the specific energy is constant up to the 9th digit after the decimal point which shows that the accuracy of the integration is quite good. A higher accuracy can be achieved by selecting smaller values for rtol and atol .","title":"Kepler Orbit"},{"location":"examples/kepler_orbit.html#kepler-orbit","text":"","title":"Kepler Orbit"},{"location":"examples/kepler_orbit.html#problem-definition","text":"For this example, we consider a spacecraft on a Kepler orbit about the Earth and we want to predict its future trajectory. We assume that the elliptical orbit is described by the following classical orbital elements: Semi-major axis: a = 20,000 km Eccentricity: e = 0.7 Inclination: i = 35\u00b0 Right ascension of the ascending node: \u03a9 = 100\u00b0 Argument of perigee: \u03c9 = 65\u00b0 True anomaly: \u03bd = 30\u00b0 The period of an orbit is given by P = 2\\pi\\sqrt{\\frac{a^3}{\\mu}} where \\mu is the standard gravitational parameter. For the Earth, \u03bc = 398600.4354 km 3 /s 2 and thus P = 2.8149 \\cdot 10 4 s = 7.82 hours. The initial position of the spacecraft expressed in Cartesian coordinates is \\mathbf{r} = \\begin{bmatrix} -5007.2484 & -1444.9181 & 3628.5346 \\end{bmatrix} \\quad km and velocity \\mathbf{v} = \\begin{bmatrix} 0.7177 & -10.2241 & 0.7482 \\end{bmatrix} \\quad km/s such that the initial state vector is \\mathbf{s} = [\\mathbf{r}, \\mathbf{v}]^{T} . The equations of motion describing the evolution of the system are \\ddot{\\mathbf{r}} = -\\frac{\\mu}{r^3}\\mathbf{r} where r is the magnitude of \\mathbf{r} . Since the crate handles first order ODEs, this system must be transformed into a state space form representation. An appropriate change of variables y_1 = r_1 , y_2 = r_2 , y_3 = r_3 , y_4 = \\dot{r}_1 = v_1 , y_5 = \\dot{r}_2 = v_2 , y_6 = \\dot{r}_3 = v_3 yields \\dot{y}_1 = y_4 \\dot{y}_2 = y_5 \\dot{y}_3 = y_6 \\dot{y}_4 = -\\frac{\\mu}{r^3}y_1 \\dot{y}_5 = -\\frac{\\mu}{r^3}y_2 \\dot{y}_6 = -\\frac{\\mu}{r^3}y_3 which can be numerically integrated.","title":"Problem Definition"},{"location":"examples/kepler_orbit.html#implementation","text":"The problem is solved using the Dopri5 method. We first import the crate and bring the required modules into scope: use ode_solvers::dopri5::*; use ode_solvers::*; Next, two type aliases are defined: type State = Vector6; type Time = f64; and a structure containing the parameters of the model is created: struct KeplerOrbit { mu: f64, } The System trait is then implemented for KeplerOrbit : impl ode_solvers::System for KeplerOrbit { // Equations of motion of the system fn system(&self, _t: Time, y: &State, dy: &mut State) { let r = (y[0] * y[0] + y[1] * y[1] + y[2] * y[2]).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = - self.mu * y[0] / r.powi(3); dy[4] = - self.mu * y[1] / r.powi(3); dy[5] = - self.mu * y[2] / r.powi(3); } } We now implement the main function. The system of ODEs is first created with the corresponding \\mu parameter. The orbital period is then computed and the initial state is initialized. The integration will be carried out using the 5th order Dormand-Prince method. Hence, a Dopri5 structure is created with the system structure, the initial time, the final time, the time increment between two outputs, the initial state vector, the relative tolerance, and the absolute tolerance as arguments. The integration is launched by calling the method integrate() and the result is stored in res . Finally, a check on the outcome of the integration is performed and the outputs can be extracted. fn main() { let system = KeplerOrbit {mu: 398600.435436} let a: f64 = 20000.0; let period = 2.0 * PI * (a.powi(3) / system.mu).sqrt(); let y0 = State::new(-5007.248417988539, -1444.918140151374, 3628.534606178356, 0.717716656891, -10.224093784269, 0.748229399696); let mut stepper = Dopri5::new(system, 0.0, 5.0 * period, 60.0, y0, 1.0e-10, 1.0e-10); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the output... // let path = Path::new(\"./outputs/kepler_orbit_dopri5.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } }","title":"Implementation"},{"location":"examples/kepler_orbit.html#results","text":"The following figure shows the trajectory resulting from the integration. As expected, the orbit is closed. Some information on the integration process is provided in the stats variable: Number of function evaluations: 62835 Number of accepted steps: 10467 Number of rejected steps: 5 In order to assess the accuracy of the integration, we look at the specific orbital energy which is a constant of the motion. This specific energy is given by \\varepsilon = \\frac{v^2}{2} - \\frac{\\mu}{r} where v is the amplitude of the velocity and r is the distance from the Earth's center to the spacecraft. In theory, this value is constant as long as no external force acts on the system (which is the case here). From the intial conditions, the specific energy of this orbit is \\varepsilon_0 = -9.9650 \\; km^2/s^2 By computing \\varepsilon for each time step and plotting the difference \\varepsilon_i - \\varepsilon_0 , we obtain the following time evolution of the error: As can be seen on this figure, the specific energy is constant up to the 9th digit after the decimal point which shows that the accuracy of the integration is quite good. A higher accuracy can be achieved by selecting smaller values for rtol and atol .","title":"Results"},{"location":"examples/three_body_problem.html","text":"Three-Body Problem Problem Definition In this problem, we consider the circular restricted three body problem in which a spacecraft evolves under the influence of the Earth and the Moon. The second order nondimensional equations of motion of this problem are \\ddot{x}-2\\dot{y}-x = -\\frac{(1-\\mu)(x+\\mu)}{d^3}-\\frac{\\mu(x-1+\\mu)}{r^3} \\ddot{y}+2\\dot{x}-y = -\\frac{(1-\\mu)y}{d^3}-\\frac{\\mu y}{r^3} \\ddot{z} = -\\frac{(1-\\mu)z}{d^3}-\\frac{\\mu z}{r^3} where d=\\left[(x+\\mu)^2+y^2+z^2\\right]^{\\frac{1}{2}} r=\\left[(x-1+\\mu)^2+y^2+z^2\\right] \\mu is a parameter of the system. For the Earth-Moon system, \\mu=0.0123 . Rewriting this system in state space form yields \\dot{x}_1=x_4 \\dot{x}_2=x_5 \\dot{x}_3=x_6 \\dot{x}_4=x_1+2x_5-\\frac{(1-\\mu)(x_1+\\mu)}{d^3}-\\frac{\\mu(x_1-1+\\mu)}{r^3} \\dot{x}_5=-2x_4+x_2-\\frac{(1-\\mu)x_2}{d^3}-\\frac{\\mu x_2}{r^3} \\dot{x}_5=-\\frac{(1-\\mu)x_3}{d^3}-\\frac{\\mu x_3}{r^3} Implementation The problem is solved using Dop853. The crate is first imported and the modules are brought into scope: use ode_solvers::dop853::*; use ode_solvers::*; Type aliases are then defined for the state vector and the time: type State = Vector6; type Time = f64; Next, we define the problem specific structure which will contain a single parameter: struct ThreeBodyProblem { mu: f64, } for which we implement the System trait: impl ode_solvers::System for ThreeBodyProblem { fn system(&self, _t: Time, y: &State, dy: &mut State) { let d = ((y[0] + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); let r = ((y[0] - 1.0 + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = y[0] + 2.0 * y[4] - (1.0 - self.mu) * (y[0] + self.mu) / d.powi(3) - self.mu * (y[0] - 1.0 + self.mu) / r.powi(3); dy[4] = -2.0 * y[3] + y[1] - (1.0 - self.mu) * y[1] / d.powi(3) - self.mu * y[1] / r.powi(3); dy[5] = -(1.0 - self.mu) * y[2] / d.powi(3) - self.mu * y[2] / r.powi(3); } } The system is created together with the intial state vector of the spacecraft in the main function. A Dop853 structure is created and the integrate() function is called. The result is handled with a match statement. fn main() { let system = ThreeBodyProblem {mu: 0.012300118882173}; let y0 = State::new(-0.271, -0.42, 0.0, 0.3, -1.0, 0.0); let mut stepper = Dop853::new(system, 0.0, 150.0, 0.001, y0, 1.0e-14, 1.0e-14); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the results... // let path = Path::new(\"./outputs/three_body_dop853.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } } Results The trajectory of the spacecraft for a duration of 150 (dimensionless time) is shown on the following figure: On that figure, the Earth is located at (0, 0) and the Moon is at (0.9877, 0). We see that the spacecraft is first on an orbit about the Earth before moving to the vicinity of the Moon. The accuracy of the integration can be assessed by looking at the Jacobi constant. The Jacobi constant is given by C=\\frac{2(1-\\mu)}{d} + \\frac{2\\mu}{r} + x^2 +y^2 - v^2 and is an integral of the motion (i.e. is constant as long as no external force is applied). Based on the initial conditions, the Jacobi constant yields C_0=3.1830 By computing the Jacobi constant at each time step and plotting the difference C-C_0 , we get the figure below: As can be seen on that figure, the Jacobi constant is constant up to the 11th digit after the decimal point. This accuracy is quite reasonable for that kind of problem but if a higher accuracy is required, the relative tolerance and/or the absolute tolerance could be reduced.","title":"Three-body Problem"},{"location":"examples/three_body_problem.html#three-body-problem","text":"","title":"Three-Body Problem"},{"location":"examples/three_body_problem.html#problem-definition","text":"In this problem, we consider the circular restricted three body problem in which a spacecraft evolves under the influence of the Earth and the Moon. The second order nondimensional equations of motion of this problem are \\ddot{x}-2\\dot{y}-x = -\\frac{(1-\\mu)(x+\\mu)}{d^3}-\\frac{\\mu(x-1+\\mu)}{r^3} \\ddot{y}+2\\dot{x}-y = -\\frac{(1-\\mu)y}{d^3}-\\frac{\\mu y}{r^3} \\ddot{z} = -\\frac{(1-\\mu)z}{d^3}-\\frac{\\mu z}{r^3} where d=\\left[(x+\\mu)^2+y^2+z^2\\right]^{\\frac{1}{2}} r=\\left[(x-1+\\mu)^2+y^2+z^2\\right] \\mu is a parameter of the system. For the Earth-Moon system, \\mu=0.0123 . Rewriting this system in state space form yields \\dot{x}_1=x_4 \\dot{x}_2=x_5 \\dot{x}_3=x_6 \\dot{x}_4=x_1+2x_5-\\frac{(1-\\mu)(x_1+\\mu)}{d^3}-\\frac{\\mu(x_1-1+\\mu)}{r^3} \\dot{x}_5=-2x_4+x_2-\\frac{(1-\\mu)x_2}{d^3}-\\frac{\\mu x_2}{r^3} \\dot{x}_5=-\\frac{(1-\\mu)x_3}{d^3}-\\frac{\\mu x_3}{r^3}","title":"Problem Definition"},{"location":"examples/three_body_problem.html#implementation","text":"The problem is solved using Dop853. The crate is first imported and the modules are brought into scope: use ode_solvers::dop853::*; use ode_solvers::*; Type aliases are then defined for the state vector and the time: type State = Vector6; type Time = f64; Next, we define the problem specific structure which will contain a single parameter: struct ThreeBodyProblem { mu: f64, } for which we implement the System trait: impl ode_solvers::System for ThreeBodyProblem { fn system(&self, _t: Time, y: &State, dy: &mut State) { let d = ((y[0] + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); let r = ((y[0] - 1.0 + self.mu).powi(2) + y[1].powi(2) + y[2].powi(2)).sqrt(); dy[0] = y[3]; dy[1] = y[4]; dy[2] = y[5]; dy[3] = y[0] + 2.0 * y[4] - (1.0 - self.mu) * (y[0] + self.mu) / d.powi(3) - self.mu * (y[0] - 1.0 + self.mu) / r.powi(3); dy[4] = -2.0 * y[3] + y[1] - (1.0 - self.mu) * y[1] / d.powi(3) - self.mu * y[1] / r.powi(3); dy[5] = -(1.0 - self.mu) * y[2] / d.powi(3) - self.mu * y[2] / r.powi(3); } } The system is created together with the intial state vector of the spacecraft in the main function. A Dop853 structure is created and the integrate() function is called. The result is handled with a match statement. fn main() { let system = ThreeBodyProblem {mu: 0.012300118882173}; let y0 = State::new(-0.271, -0.42, 0.0, 0.3, -1.0, 0.0); let mut stepper = Dop853::new(system, 0.0, 150.0, 0.001, y0, 1.0e-14, 1.0e-14); let res = stepper.integrate(); // Handle result match res { Ok(stats) => { stats.print(); // Do something with the results... // let path = Path::new(\"./outputs/three_body_dop853.dat\"); // save(stepper.x_out(), stepper.y_out(), path); // println!(\"Results saved in: {:?}\", path); }, Err(_) => println!(\"An error occured.\"), } }","title":"Implementation"},{"location":"examples/three_body_problem.html#results","text":"The trajectory of the spacecraft for a duration of 150 (dimensionless time) is shown on the following figure: On that figure, the Earth is located at (0, 0) and the Moon is at (0.9877, 0). We see that the spacecraft is first on an orbit about the Earth before moving to the vicinity of the Moon. The accuracy of the integration can be assessed by looking at the Jacobi constant. The Jacobi constant is given by C=\\frac{2(1-\\mu)}{d} + \\frac{2\\mu}{r} + x^2 +y^2 - v^2 and is an integral of the motion (i.e. is constant as long as no external force is applied). Based on the initial conditions, the Jacobi constant yields C_0=3.1830 By computing the Jacobi constant at each time step and plotting the difference C-C_0 , we get the figure below: As can be seen on that figure, the Jacobi constant is constant up to the 11th digit after the decimal point. This accuracy is quite reasonable for that kind of problem but if a higher accuracy is required, the relative tolerance and/or the absolute tolerance could be reduced.","title":"Results"}]}
\ No newline at end of file
diff --git a/sitemap.xml.gz b/sitemap.xml.gz
index d7c63f5778b01e90224e7a43835d5ea01229a69f..d7259fb1d5dc4215272e938f02a16a006d31d956 100644
GIT binary patch
delta 14
VcmeBR>R@7(@8;mRx{;Nc5da$+17iRH

delta 14
VcmeBR>R@7(@8;kz+Q`bx2mlxB0>A(O