From 71080f8c139ee5b9463cabe65166e5c621de408b Mon Sep 17 00:00:00 2001 From: John Long Date: Wed, 10 Jan 2024 13:14:10 -0500 Subject: [PATCH] Made README more detailed, added docstrings --- README.md | 62 ++++++++++++++++- src/interface.jl | 81 ++++++++++++++++++++++ src/types.jl | 170 ++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 310 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index fd1b9c3..17978bf 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ [![codecov](https://codecov.io/github/QuEraComputing/DormandPrince.jl/graph/badge.svg?token=qYZ4V7m0JY)](https://codecov.io/github/QuEraComputing/DormandPrince.jl) -Julia-native implementation of the Dormand-Prince 5th order solver +Julia-native implementation of the Dormand-Prince 5th and 8th order solvers ## Installation @@ -22,6 +22,66 @@ DormandPrince is a   pkg> add DormandPrince ``` +## Usage + +### Single End-Time + +```julia +julia> using DormandPrince + +# define your ODE, in this case, dy/dx = 0.85y +julia> function fcn(x, y, f) + f[1] = 0.85y[1] + end + +# Create a solver object. We will use the 5th order solver and +# start integrating at t = 0.0 with initial value of 19.0 +julia> solver = DP5Solver(fcn, 0.0, [19.0]) + +# Begin integration up to t = 2.1 +julia> integrate!(solver, 2.1) + +# get_current_state will return the answer +julia> get_current_state(solver) +``` + +Both the `DP5Solver` and `DP8Solver`'s are stateful allowing for memory-efficient integration to future end times from the last integrated end point (e.g. if you chose t = 1.0 as your endpoint you can call `integrate!` again with t=2.0 and it will "carry forward" the work starting from t = 1.0 instead of requiring you to set things up all over again). + +### Multiple End-Times + +```julia +julia> using DormandPrince + +# Define ODE +julia> function fcn(x, y, f) + f[1] = 0.85y[1] + end + +# Define times of interest to analyze/perform actions on the solution +julia> times = [1.0, 1.1, 1.9, 2.4] + +# Create the solver object, with integration starting from t = 0.0 and initial value of 19.0 +julia> solver = DP5Solver(fcn, 0.0, [19.0]) + +# Empty array to store intermediate values +julia> intermediate_values = [] + +# Use callback to store intermediate values. The solver object will also be mutated to store the solution +# found at the last time point. +julia> integrate!(solver, times) do time, val + push!(intermediate_values, val[]) + end + +``` + +You may also create a `SolverIterator` that can then be iterated over. Note that this will require a fresh solver + +```julia +for (time, value) in SolverIterator(solver, times) + println("Time: ", time, " ", "Value: ", value[]) +end +``` + ## License Apache License 2.0 diff --git a/src/interface.jl b/src/interface.jl index 229ba0c..d335926 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -1,4 +1,28 @@ +""" + struct SolverIterator{T <: Real} + SolverIterator(solver, times) + +Given a solver and a vector of times, +this iterator will return the state of the solver at each time. +The solver will be mutated to hold the solution from the last time given. + +# Examples + +```julia +julia> solver = DP5Solver(fcn, 0.0, [0.0]) + +julia> times = [1.0, 2.0, 3.0] + +julia> intermediate_vals = [] + +julia> for (time, value) in SolverIterator(solver, times) + push!(intermediate_values, value[]) + end + +``` + +""" struct SolverIterator{T <: Real} solver::AbstractDPSolver{T} times::AbstractVector{T} @@ -32,12 +56,69 @@ end get_current_state(::AbstractDPSolver) = error("not implemented") integrate_core!(::AbstractDPSolver{T}, ::T) where T = error("not implemented") +""" + integrate!(solver, time) + +Integrate the ODE problem that is part of solver to the end time `time`. +`solver` can be a `DP5Solver` or a `DP8Solver` type. + +The `solver` is mutated to hold the solution and `solver` state at the time `time`, allowing for +efficient integration to future end times of interest through subsequent calls to `integrate!` with +later times. + +# Examples + +```julia +julia> function fcn(x, y, f) + f[1] = 0.85y[1] + end + +julia> solver = DP5Solver(fcn, 0.0, [19.0]) + +# Integrate from initial time of 0.0 to 1.0 +julia> integrate!(solver, 1.0) + +# integrate from the last time of 1.0 to 2.0 +julia> integrate!(solver, 2.0) + +# Solver object now holds solution and state at 2.0 +julia> get_current_state(solver) +``` +""" function integrate!(solver::AbstractDPSolver{T}, time::T) where T <: Real report = integrate_core!(solver, time) if report.idid != COMPUTATION_SUCCESSFUL error("integration failed at time $time with report $report") end end + +""" + integrate!(callback, solver, times) + +Integrate the ODE problem that is part of `solver` to the end times `times` and apply the callback on +the time and solution at each time. +`solver` can be a `DP5Solver` or a `DP8Solver` type. + +At the end of all the times the solver holds the solution and solver state at the last time in `times`. + +# Examples + +```julia +julia> function fcn(x, y, f) + f[1] = 0.85y[1] + end + +julia> times = [1.0, 1.1, 1.9, 2.4] + +julia> solver = DP5Solver(fcn, 0.0, [19.0]) + +julia> intermediate_values = [] + +julia> integrate!(solver, times) do time, val + push!(intermediate_values, val[]) + end +``` +""" function integrate!(callback, solver::AbstractDPSolver{T}, times::AbstractVector{T}; sort_times::Bool = true) where {T <: Real} times = sort_times ? sort(collect(times)) : times diff --git a/src/types.jl b/src/types.jl index 6eda7dc..3f3ef2b 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,6 +1,19 @@ abstract type AbstractDPSolver{T <: Real, StateType <: AbstractVector, F} end +""" + @enum Idid + +Enum used to represent status of the integration process in the `Report` type. + +# Values +- `COMPUTATION_SUCCESSFUL`: Integration completed successfully. +- `INPUT_NOT_CONSISTENT`: The options given to the solver violate acceptable ranges (see the +`checks` field of the `Report` type for more information). +- `LARGER_NMAX_NEEDED`: The maximum number of allowed steps is too small. +- `STEP_SIZE_BECOMES_TOO_SMALL`: The step size becomes too small. +- `STEP_SIZE_BECOMES_NAN`: The step size becomes NaN. +""" @enum Idid begin COMPUTATION_SUCCESSFUL = 1 INPUT_NOT_CONSISTENT = -1 # use for check failures in the beginning of core_integrator call @@ -9,6 +22,18 @@ abstract type AbstractDPSolver{T <: Real, StateType <: AbstractVector, F} end STEP_SIZE_BECOMES_NAN = -4 end +""" + @enum Checks + +Enum used to represent status of checks on user options to the solver in the `Report` type. + +# Values +- `INPUT_CHECKS_SUCCESSFUL`: All checks on user options passed. +- `MAX_ALLOWED_STEPS_NEGATIVE`: The maximum allowed steps value is negative. +- `UNSUPPORTED_UROUND`: The uround value is either too small (<= 1e-35) or too large (>= 1.0). +- `CURIOUS_BETA`: The Beta value for stabilized step control is greater than 0.2. +- `CURIOUS_SAFETY_FACTOR`: The safety factor is either too large (>= 1.0) or too small (<= 1e-4). +""" @enum Checks begin INPUT_CHECKS_SUCCESSFUL MAX_ALLOWED_STEPS_NEGATIVE @@ -17,6 +42,22 @@ end CURIOUS_SAFETY_FACTOR end +""" + struct Report{T <: Real} + +Contains data on the integration process including after integration has completed and if an error was detected within +options provided for integration. + +# Fields +- `x_final`: The final time point integration reached. Should be equal to the time point provided by the user with successful usage. +- `checks`: Status of checks performed on the options provided by the user, represented by a `Checks` element. Should be `INPUT_CHECKS_SUCCESSFUL` with successful usage. +- `idid`: Status of the integration process, represented by an `Idid` element. Should be `COMPUTATION_SUCCESSFUL` with successful usage. +- `num_func_evals`: Number of function evaluations performed during integration. +- `num_computed_steps`: Number of computed steps during integration. +- `num_accpeted_steps`: Number of accepted steps during integration. +- `num_rejected_steps`: Number of rejected steps during integration. +- +""" struct Report{T <: Real} x_final::T checks::Checks @@ -28,6 +69,27 @@ struct Report{T <: Real} num_rejected_steps::Int end +""" + struct Options{T <: Real} + +Holds the options used in the integration process by a solver object. + +# Fields +- `atol`: Absolute Tolerance. Default is 1e-10 +- `rtol`: Relative Tolerance. Default is 1e-10 +- `uround`: Rounding unit, default is eps(T) where T <: Real. +- `safety_factor`: Safety factor in step size prediction, default is 0.9. +- Step Size Selection parameters, with the new step size being subject to the constraint: `step_size_selection_one <= new_step/old_step <= step_size_selection_two` + - `step_size_selection_one`: Defaut is 0.2 + - `step_size_selection_two`: Default is 10.0 +- `beta`: Beta for stabilized step size control. Large values of beta (<= 0.1) make step size control more stable. +- `maximal_step_size`: The largest the step size can be. Default is 0.0, which internally translates to `xend - x`. +- `initial_step_size`: Initial step size, default is 0.0. An initial guess is computed internally. +- `maximum_allowed_steps`: Maximum number of allowed steps, default is 100000. +- `print_error_messages`: Whether to print error messages, default is true. +- `stiffness_test_activation_step`: After integer multiples of this number of steps, perform stiffness detection. Default is 1000. + +""" @option struct Options{T <: Real} # originally in work[1] - work[7] uround::T = eps(T) @@ -78,8 +140,19 @@ end end -# should " core_integrator" take in DP5Solver or should DP5Solver have some associated method -# attached to it? +""" + struct DP5Solver + +A 5th Order Dormand-Prince solver object that contains: +- the ODE problem of interest +- The solution to the ODE at the last integrated-to time point +- intermediate arrays used in the Runge-Kutta method +- constants and variables used by the solver +- user-provided options for the solver + +The storage of such intermediate information allows for efficient integration from a previously integrated-to +time point and a future time point. +""" struct DP5Solver{T, StateType , F, OptionsType, ConstsType, VarsType} <: AbstractDPSolver{T, StateType, F} f::F y::StateType @@ -120,6 +193,41 @@ struct DP5Solver{T, StateType , F, OptionsType, ConstsType, VarsType} <: Abstrac end end +""" + DP5Solver(f, x, y; kw...) + +Create a 5th Order Dormand-Prince solver object to solve the ODE problem `y' = f(x, y)`. + +# Examples + +```julia +julia> fcn(x, y, f) + f[1] = 0.85y[1] + end + +julia> solver = DP5Solver(fcn, 0.0, [19.0]; atol = 1e-10, rtol = 1e-10) +``` + +# Arguments +- `f`: The function representing the ODE, should be in the form `f(x, y, f)`. +- `x`: The starting time point of the ODE problem. +- `y`: The initial value of the ODE in vector form + +# Keyword Arguments +- `atol`: Absolute Tolerance. Default is 1e-10 +- `rtol`: Relative Tolerance. Default is 1e-10 +- `uround`: Rounding unit, default is eps(T) where T <: Real. +- `safety_factor`: Safety factor in step size prediction, default is 0.9. +- Step Size Selection parameters, with the new step size being subject to the constraint: `step_size_selection_one <= new_step/old_step <= step_size_selection_two` + - `step_size_selection_one`: Defaut is 0.2 + - `step_size_selection_two`: Default is 10.0 +- `beta`: Beta for stabilized step size control. Large values of beta (<= 0.1) make step size control more stable. +- `maximal_step_size`: The largest the step size can be. Default is 0.0, which internally translates to `xend - x`. +- `initial_step_size`: Initial step size, default is 0.0. An initial guess is computed internally. +- `maximum_allowed_steps`: Maximum number of allowed steps, default is 100000. +- `print_error_messages`: Whether to print error messages, default is true. +- `stiffness_test_activation_step`: After integer multiples of this number of steps, perform stiffness detection. Default is 1000. +""" function DP5Solver( f, x::Real, @@ -138,6 +246,19 @@ end # should " core_integrator" take in DP5Solver or should DP5Solver have some associated method # attached to it? +""" + struct DP8Solver + +An 8th Order Dormand-Prince solver object that contains: +- the ODE problem of interest +- The solution to the ODE at the last integrated-to time point +- intermediate arrays used in the Runge-Kutta method +- constants and variables used by the solver +- user-provided options for the solver + +The storage of such intermediate information allows for efficient integration from a previously integrated-to +time point and a future time point. +""" struct DP8Solver{T, StateType ,F, OptionsType, ConstsType, VarsType} <: AbstractDPSolver{T, StateType, F} f::F y::StateType @@ -195,6 +316,41 @@ struct DP8Solver{T, StateType ,F, OptionsType, ConstsType, VarsType} <: Abstract end end +""" + DP8Solver(f, x, y; kw...) + +Create an 8th Order Dormand-Prince solver object to solve the ODE problem `y' = f(x, y)`. + +# Examples + +```julia +julia> fcn(x, y, f) + f[1] = 0.85y[1] + end + +julia> solver = DP8Solver(fcn, 0.0, [19.0]; atol = 1e-10, rtol = 1e-10) +``` + +# Arguments +- `f`: The function representing the ODE, should be in the form `f(x, y, f)`. +- `x`: The starting time point of the ODE problem. +- `y`: The initial value of the ODE in vector form + +# Keyword Arguments +- `atol`: Absolute Tolerance. Default is 1e-10 +- `rtol`: Relative Tolerance. Default is 1e-10 +- `uround`: Rounding unit, default is eps(T) where T <: Real. +- `safety_factor`: Safety factor in step size prediction, default is 0.9. +- Step Size Selection parameters, with the new step size being subject to the constraint: `step_size_selection_one <= new_step/old_step <= step_size_selection_two` + - `step_size_selection_one`: Defaut is 0.2 + - `step_size_selection_two`: Default is 10.0 +- `beta`: Beta for stabilized step size control. Large values of beta (<= 0.1) make step size control more stable. +- `maximal_step_size`: The largest the step size can be. Default is 0.0, which internally translates to `xend - x`. +- `initial_step_size`: Initial step size, default is 0.0. An initial guess is computed internally. +- `maximum_allowed_steps`: Maximum number of allowed steps, default is 100000. +- `print_error_messages`: Whether to print error messages, default is true. +- `stiffness_test_activation_step`: After integer multiples of this number of steps, perform stiffness detection. Default is 1000. +""" function DP8Solver( f, x::Real, @@ -213,5 +369,15 @@ function DP8Solver( DP8Solver(f, x, y, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, y1;kw...) end +""" + get_current_state(solver::DP5Solver) + +Gets the current state/solution held by the `DP5Solver` object. +""" get_current_state(solver::DP5Solver) = solver.y +""" + get_current_state(solver::DP8Solver) + +Gets the current state/solution held by the `DP8Solver` object. +""" get_current_state(solver::DP8Solver) = solver.y \ No newline at end of file