diff --git a/previews/PR194/.documenter-siteinfo.json b/previews/PR194/.documenter-siteinfo.json index b0724c9ea..043315239 100644 --- a/previews/PR194/.documenter-siteinfo.json +++ b/previews/PR194/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-07-12T10:25:29","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-07-12T10:26:52","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/previews/PR194/api/index.html b/previews/PR194/api/index.html index 6742e823d..d18aac2b5 100644 --- a/previews/PR194/api/index.html +++ b/previews/PR194/api/index.html @@ -16,14 +16,14 @@ tv, xv = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100 ) -tv, xv, psol = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100, Val(true) )source
TaylorIntegration.lyap_taylorinteg — Function
lyap_taylorinteg(f!, q0, t0, tmax, order, abstol[, jacobianfunc!=nothing];
+tv, xv, psol = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100, Val(true) )
source
TaylorIntegration.lyap_taylorinteg — Function
lyap_taylorinteg(f!, q0, t0, tmax, order, abstol[, jacobianfunc!=nothing];
     maxsteps::Int=500, parse_eqs::Bool=true)
 lyap_taylorinteg(f!, q0, trange, order, abstol[, jacobianfunc!=nothing];
-    maxsteps::Int=500, parse_eqs::Bool=true)

Similar to taylorinteg for the calculation of the Lyapunov spectrum. Note that the number of TaylorN variables should be set previously by the user (e.g., by means of TaylorSeries.set_variables) and should be equal to the length of the vector of initial conditions q0. Otherwise, whenever length(q0) != TaylorSeries.get_numvars(), then lyap_taylorinteg throws an AssertionError. Optionally, the user may provide a Jacobian function jacobianfunc! to evaluate the current value of the Jacobian. Otherwise, the current value of the Jacobian is computed via automatic differentiation using TaylorSeries.jl.

source
TaylorIntegration.@taylorize — Macro

@taylorize expr

This macro evals the function given by expr and defines a new method of jetcoeffs! which is specialized on that function. Integrating via taylorinteg of lyap_taylorinteg after using the macro yields better performance.

See the documentation for more details and limitations.

Warning

This macro is on an experimental stage; check the integration results carefully.

source

Internal

TaylorIntegration.BookKeeping — Type
BookKeeping

Mutable struct that contains all the bookkeeping vectors/dictionaries used within _make_parsed_jetcoeffs: - d_indx : Dictionary mapping new variables (symbols) to old (perhaps indexed) symbols - d_assign : Dictionary with the numeric assignments (that are substituted) - d_decl : Dictionary declared arrays - v_newvars : Symbols of auxiliary indexed vars - v_arraydecl: Symbols which are explicitly declared as Array or Vector - v_array1 : Symbols which are explicitly declared as Array{Taylor1{T},1} - v_array2 : Symbols which are explicitly declared as Array{Taylor1{T},2} - v_array3 : Symbols which are explicitly declared as Array{Taylor1{T},3} - v_array4 : Symbols which are explicitly declared as Array{Taylor1{T},4} - v_preamb : Symbols or Expr used in the preamble (declarations, etc) - retvar : Guessed returned variable, which defines the LHS of the ODEs

source
TaylorIntegration.RetAlloc — Type
RetAlloc{Taylor1{T}}

Struct related to the returned variables that are pre-allocated when @taylorize is used. - v0 : Vector{Taylor1{T}} - v1 : Vector{Vector{Taylor1{T}}}

source
TaylorIntegration.__jetcoeffs! — Method
__jetcoeffs!(::Val{false}, f, t, x, params)
+    maxsteps::Int=500, parse_eqs::Bool=true)

Similar to taylorinteg for the calculation of the Lyapunov spectrum. Note that the number of TaylorN variables should be set previously by the user (e.g., by means of TaylorSeries.set_variables) and should be equal to the length of the vector of initial conditions q0. Otherwise, whenever length(q0) != TaylorSeries.get_numvars(), then lyap_taylorinteg throws an AssertionError. Optionally, the user may provide a Jacobian function jacobianfunc! to evaluate the current value of the Jacobian. Otherwise, the current value of the Jacobian is computed via automatic differentiation using TaylorSeries.jl.

source
TaylorIntegration.@taylorize — Macro

@taylorize expr

This macro evals the function given by expr and defines a new method of jetcoeffs! which is specialized on that function. Integrating via taylorinteg of lyap_taylorinteg after using the macro yields better performance.

See the documentation for more details and limitations.

Warning

This macro is on an experimental stage; check the integration results carefully.

source

Internal

TaylorIntegration.BookKeeping — Type
BookKeeping

Mutable struct that contains all the bookkeeping vectors/dictionaries used within _make_parsed_jetcoeffs: - d_indx : Dictionary mapping new variables (symbols) to old (perhaps indexed) symbols - d_assign : Dictionary with the numeric assignments (that are substituted) - d_decl : Dictionary declared arrays - v_newvars : Symbols of auxiliary indexed vars - v_arraydecl: Symbols which are explicitly declared as Array or Vector - v_array1 : Symbols which are explicitly declared as Array{Taylor1{T},1} - v_array2 : Symbols which are explicitly declared as Array{Taylor1{T},2} - v_array3 : Symbols which are explicitly declared as Array{Taylor1{T},3} - v_array4 : Symbols which are explicitly declared as Array{Taylor1{T},4} - v_preamb : Symbols or Expr used in the preamble (declarations, etc) - retvar : Guessed returned variable, which defines the LHS of the ODEs

source
TaylorIntegration.RetAlloc — Type
RetAlloc{Taylor1{T}}

Struct related to the returned variables that are pre-allocated when @taylorize is used. - v0 : Vector{Taylor1{T}} - v1 : Vector{Vector{Taylor1{T}}}

source
TaylorIntegration.__jetcoeffs! — Method
__jetcoeffs!(::Val{false}, f, t, x, params)
 __jetcoeffs!(::Val{true}, f, t, x, params, rv)
 __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params)
-__jetcoeffs!(::Val{true}, f, t, x, dx, params, rv)

Chooses a method of jetcoeffs! (hard-coded) or the generated by @taylorize) depending on Val{bool} (bool::Bool).

source
TaylorIntegration._allocated_defs! — Method

_allocated_defs!(new_jetcoeffs, bkkeep)

Add allocated variable definitions to new_jetcoeffs, to make it more human readable.

source
TaylorIntegration._capture_fn_args_body! — Function

_capture_fn_args_body!(ex, dd::Dict{Symbol, Any})

Captures the name of a function, arguments, body and other properties, returning them as the values of the dictionary dd, which is updated in place.

source
TaylorIntegration._defs_allocs! — Function

_defs_allocs!(preamble, fnargs, bkkeep, [inloop=false, ex_aux::Expr(:block,)])

Returns a vector with expressions defining the auxiliary variables in the preamble, and the declaration of the arrays. This function may modify bkkeep.d_indx if new variables are introduced. bkkeep.v_preamb is for bookkeeping the introduced variables.

source
TaylorIntegration._determine_parsing! — Method
_determine_parsing!(parse_eqs::Bool, f, t, x, params)
-_determine_parsing!(parse_eqs::Bool, f, t, x, dx, params)

Check if the parsed method of jetcoeffs! exists and check it runs without error.

source
TaylorIntegration._extract_parts — Method

_extract_parts(ex::Expr)

Returns the function name, the function arguments, and the body of a function passed as an Expr. The function may be provided as a one-line function, or in the long form (anonymous functions do not work).

source
TaylorIntegration._make_parsed_jetcoeffs — Method

_make_parsed_jetcoeffs( ex )

This function constructs the expressions of two new methods, the first equivalent to the differential equations (jetcoeffs!), which exploits the mutating functions of TaylorSeries.jl, and the second one (allocatejetcoeffs) preallocates any auxiliary Taylor1 or Vector{Taylor1{T}} needed.

source
TaylorIntegration._newfnbody! — Method

_newfnbody!(fnbody, fnargs, bkkeep)

Returns a new (modified) body of the function, a priori unfolding the expression graph (AST) as unary and binary calls, and updates the bookkeeping structure bkkeep.

source
TaylorIntegration._newhead — Method

_newhead(fn, fnargs)

Creates the head of the new method of jetcoeffs! and _allocate_jetcoeffs. fn is the name of the passed function and fnargs is a vector with its arguments defning the function (which are either three or four).

source
TaylorIntegration._parse_newfnbody! — Method

parsenewfnbody!(ex::Expr, preex::Expr, prealloc::Expr, bkkeep::BookKeeping, inloop::Bool)

Parses ex (the new body of the function) replacing the expressions to use the mutating functions of TaylorSeries, and building the preamble preex and prealloc expressions. This is done by traversing recursively (again) the args of ex, updating the bookkeeping struct bkkeep, in particular the fieldnames v_newvars and d_assign.

source
TaylorIntegration._preamble_body — Method

_preamble_body(fnbody, fnargs)

Returns expressions for the preamble, the declaration of arrays, the body and the bookkeeping struct, which will be used to build the new functions. fnbody is the expression with the body of the original function (already adapted), fnargs is a vector of symbols of the original diferential equations function.

source
TaylorIntegration._recursionloop — Method

_recursionloop(fnargs, bkkeep)

Build the expression for the recursion-loop.

source
TaylorIntegration._rename_indexedvars — Method

_rename_indexedvars(fnbody)

Renames the indexed variables (using Espresso.genname()) that exists in fnbody. Returns fnbody with the renamed variables and a dictionary that links the new variables to the old indexed ones.

source
TaylorIntegration._replace_expr! — Method

_replace_expr!(ex::Expr, preex::Expr, , prealloc::Expr, i::Int, aalhs, aarhs, bkkeep::BookKeeping)

Replaces the calls in ex.args[i], and updates preex and prealloc with the appropriate expressions, based on the the LHS (aalhs) and RHS (aarhs) of the base assignment. The bookkeeping struct is updated (v_newvars) within _replacecalls!. d_indx is used to bring back the indexed variables.

source
TaylorIntegration._replacecalls! — Method

_replacecalls!(bkkeep, fnold, newvar)

Replaces the symbols of unary and binary calls of the expression fnold, which defines newvar, by the mutating functions in TaylorSeries.jl. The vector bkkeep.v_vars is updated if new auxiliary variables are introduced (bookkeeping).

source
TaylorIntegration._returned_expr — Method

_returned_expr(bkkeep)

Constructs the expression to be returned by TaylorIntegration._allocate_jetcoeffs!

source
TaylorIntegration._second_stepsize — Method
_second_stepsize(x, epsilon)

Corresponds to the "second stepsize control" in Jorba and Zou (2005) paper. We use it if stepsize returns Inf.

source
TaylorIntegration._split_arraydecl! — Method

_split_arraydecl!(bkkeep)

Split bkkeep.varraydecl in the vector (bkkeep.varray1), matrix (bkkeep.v_array2), etc, to properly construct the RetAlloc variable.

source
TaylorIntegration._stepsize — Method
_stepsize(aux1, epsilon, k)

Helper function to avoid code repetition. Returns $(epsilon/aux1)^(1/k)$.

source
TaylorIntegration.findroot! — Method
findroot!(t, x, dx, g_tupl_old, g_tupl, eventorder, tvS, xvS, gvS,
+__jetcoeffs!(::Val{true}, f, t, x, dx, params, rv)

Chooses a method of jetcoeffs! (hard-coded) or the generated by @taylorize) depending on Val{bool} (bool::Bool).

source
TaylorIntegration._allocated_defs! — Method

_allocated_defs!(new_jetcoeffs, bkkeep)

Add allocated variable definitions to new_jetcoeffs, to make it more human readable.

source
TaylorIntegration._capture_fn_args_body! — Function

_capture_fn_args_body!(ex, dd::Dict{Symbol, Any})

Captures the name of a function, arguments, body and other properties, returning them as the values of the dictionary dd, which is updated in place.

source
TaylorIntegration._defs_allocs! — Function

_defs_allocs!(preamble, fnargs, bkkeep, [inloop=false, ex_aux::Expr(:block,)])

Returns a vector with expressions defining the auxiliary variables in the preamble, and the declaration of the arrays. This function may modify bkkeep.d_indx if new variables are introduced. bkkeep.v_preamb is for bookkeeping the introduced variables.

source
TaylorIntegration._determine_parsing! — Method
_determine_parsing!(parse_eqs::Bool, f, t, x, params)
+_determine_parsing!(parse_eqs::Bool, f, t, x, dx, params)

Check if the parsed method of jetcoeffs! exists and check it runs without error.

source
TaylorIntegration._extract_parts — Method

_extract_parts(ex::Expr)

Returns the function name, the function arguments, and the body of a function passed as an Expr. The function may be provided as a one-line function, or in the long form (anonymous functions do not work).

source
TaylorIntegration._make_parsed_jetcoeffs — Method

_make_parsed_jetcoeffs( ex )

This function constructs the expressions of two new methods, the first equivalent to the differential equations (jetcoeffs!), which exploits the mutating functions of TaylorSeries.jl, and the second one (allocatejetcoeffs) preallocates any auxiliary Taylor1 or Vector{Taylor1{T}} needed.

source
TaylorIntegration._newfnbody! — Method

_newfnbody!(fnbody, fnargs, bkkeep)

Returns a new (modified) body of the function, a priori unfolding the expression graph (AST) as unary and binary calls, and updates the bookkeeping structure bkkeep.

source
TaylorIntegration._newhead — Method

_newhead(fn, fnargs)

Creates the head of the new method of jetcoeffs! and _allocate_jetcoeffs. fn is the name of the passed function and fnargs is a vector with its arguments defning the function (which are either three or four).

source
TaylorIntegration._parse_newfnbody! — Method

parsenewfnbody!(ex::Expr, preex::Expr, prealloc::Expr, bkkeep::BookKeeping, inloop::Bool)

Parses ex (the new body of the function) replacing the expressions to use the mutating functions of TaylorSeries, and building the preamble preex and prealloc expressions. This is done by traversing recursively (again) the args of ex, updating the bookkeeping struct bkkeep, in particular the fieldnames v_newvars and d_assign.

source
TaylorIntegration._preamble_body — Method

_preamble_body(fnbody, fnargs)

Returns expressions for the preamble, the declaration of arrays, the body and the bookkeeping struct, which will be used to build the new functions. fnbody is the expression with the body of the original function (already adapted), fnargs is a vector of symbols of the original diferential equations function.

source
TaylorIntegration._recursionloop — Method

_recursionloop(fnargs, bkkeep)

Build the expression for the recursion-loop.

source
TaylorIntegration._rename_indexedvars — Method

_rename_indexedvars(fnbody)

Renames the indexed variables (using Espresso.genname()) that exists in fnbody. Returns fnbody with the renamed variables and a dictionary that links the new variables to the old indexed ones.

source
TaylorIntegration._replace_expr! — Method

_replace_expr!(ex::Expr, preex::Expr, , prealloc::Expr, i::Int, aalhs, aarhs, bkkeep::BookKeeping)

Replaces the calls in ex.args[i], and updates preex and prealloc with the appropriate expressions, based on the the LHS (aalhs) and RHS (aarhs) of the base assignment. The bookkeeping struct is updated (v_newvars) within _replacecalls!. d_indx is used to bring back the indexed variables.

source
TaylorIntegration._replacecalls! — Method

_replacecalls!(bkkeep, fnold, newvar)

Replaces the symbols of unary and binary calls of the expression fnold, which defines newvar, by the mutating functions in TaylorSeries.jl. The vector bkkeep.v_vars is updated if new auxiliary variables are introduced (bookkeeping).

source
TaylorIntegration._returned_expr — Method

_returned_expr(bkkeep)

Constructs the expression to be returned by TaylorIntegration._allocate_jetcoeffs!

source
TaylorIntegration._second_stepsize — Method
_second_stepsize(x, epsilon)

Corresponds to the "second stepsize control" in Jorba and Zou (2005) paper. We use it if stepsize returns Inf.

source
TaylorIntegration._split_arraydecl! — Method

_split_arraydecl!(bkkeep)

Split bkkeep.varraydecl in the vector (bkkeep.varray1), matrix (bkkeep.v_array2), etc, to properly construct the RetAlloc variable.

source
TaylorIntegration._stepsize — Method
_stepsize(aux1, epsilon, k)

Helper function to avoid code repetition. Returns $(epsilon/aux1)^(1/k)$.

source
TaylorIntegration.findroot! — Method
findroot!(t, x, dx, g_tupl_old, g_tupl, eventorder, tvS, xvS, gvS,
     t0, δt_old, x_dx, x_dx_val, g_dg, g_dg_val, nrabstol,
-    newtoniter, nevents) -> nevents

Internal root-finding subroutine, based on Newton-Raphson process. If there is a crossing, then the crossing data is stored in tvS, xvS and gvS and nevents, the number of events/crossings, is updated. Here t is a Taylor1 polynomial which represents the independent variable; x is an array of Taylor1 variables which represent the vector of dependent variables; dx is an array of Taylor1 variables which represent the LHS of the ODE; g_tupl_old is the last-before-current value returned by event function g and g_tupl is the current one; eventorder is the order of the derivative of g whose roots the user is interested in finding; tvS stores the surface-crossing instants; xvS stores the value of the solution at each of the crossings; gvS stores the values of the event function g (or its eventorder-th derivative) at each of the crossings; t0 is the current time; δt_old is the last time-step size; x_dx, x_dx_val, g_dg, g_dg_val are auxiliary variables; nrabstol is the Newton-Raphson process tolerance; newtoniter is the maximum allowed number of Newton-Raphson iteration; nevents is the current number of detected events/crossings.

source
TaylorIntegration.inbookkeeping — Method

inbookkeeping(v, bkkeep::BookKeeping)

Checks if v is declared in bkkeep, considering the d_indx, v_newvars and v_arraydecl fields.

source
TaylorIntegration.jetcoeffs! — Method
jetcoeffs!(eqsdiff!::Function, t, x, dx, xaux, params)

Mutates x in-place using the recursion relation of the derivatives obtained from the differential equations $\dot{x}=dx/dt=f(x, p, t)$.

eqsdiff! is the function defining the RHS of the ODE, x contains the Taylor1 expansion of the dependent variables and t is the independent variable, and params are the parameters appearing on the function defining the differential equation. See taylorinteg for examples and convention for eqsdiff. Note that x is of type Vector{Taylor1{U}} where U<:Number; t is of type Taylor1{T} where T<:Real. In this case, two auxiliary containers dx and xaux (both of the same type as x) are needed to avoid allocations.

Initially, x contains only the 0-th order Taylor coefficient of the current system state (the initial conditions), and jetcoeffs! computes recursively the high-order derivates back into x.

source
TaylorIntegration.jetcoeffs! — Method
jetcoeffs!(eqsdiff::Function, t, x, params)

Returns an updated x using the recursion relation of the derivatives obtained from the differential equations $\dot{x}=dx/dt=f(x, p, t)$.

eqsdiff is the function defining the RHS of the ODE, x contains the Taylor1 expansion of the dependent variable(s) and t is the independent variable, and params are the parameters appearing on the function defining the differential equation. See taylorinteg for examples and convention for eqsdiff. Note that x is of type Taylor1{U} where U<:Number; t is of type Taylor1{T} where T<:Real.

Initially, x contains only the 0-th order Taylor coefficient of the current system state (the initial conditions), and jetcoeffs! computes recursively the high-order derivates back into x.

source
TaylorIntegration.lyap_jetcoeffs! — Method
lyap_jetcoeffs!(t, x, dx, jac, varsaux)

Similar to jetcoeffs! for the calculation of the Lyapunov spectrum. Updates only the elements of x which correspond to the solution of the 1st-order variational equations $\dot{\xi}=J \cdot \xi$, where $J$ is the Jacobian matrix, i.e., the linearization of the equations of motion. jac is the Taylor expansion of $J$ wrt the independent variable, around the current initial condition. varsaux is an auxiliary array of type Array{eltype(jac),3} to avoid allocations. Calling this method assumes that jac has been computed previously using stabilitymatrix!.

source
TaylorIntegration.lyap_taylorstep! — Method
lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, t0, t1, order, abstol, _δv, varsaux, params[, jacobianfunc!])

Similar to taylorstep! for the calculation of the Lyapunov spectrum. jac is the Taylor expansion (wrt the independent variable) of the linearization of the equations of motion, i.e, the Jacobian. xaux, δx, dδx, varsaux and _δv are auxiliary vectors, and params define the parameters of the ODEs. Optionally, the user may provide a Jacobian function jacobianfunc! to compute jac. Otherwise, jac is computed via automatic differentiation using TaylorSeries.jl.

source
TaylorIntegration.nrconvergencecriterion — Method
nrconvergencecriterion(g_val, nrabstol::T, nriter::Int, newtoniter::Int) where {T<:Real}

A rudimentary convergence criterion for the Newton-Raphson root-finding process. g_val may be either a Real, Taylor1{T} or a TaylorN{T}, where T<:Real. Returns true if: 1) the absolute value of g_val, the value of the event function g evaluated at the current estimated root by the Newton-Raphson process, is less than the nrabstol tolerance; and 2) the number of iterations nriter of the Newton-Raphson process is less than the maximum allowed number of iterations, newtoniter; otherwise, returns false.

source
TaylorIntegration.stabilitymatrix! — Method
stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, _δv, params[, jacobianfunc!=nothing])

Updates the matrix jac::Matrix{Taylor1{U}} (linearized equations of motion) computed from the equations of motion (eqsdiff!), at time t at x; x is of type Vector{Taylor1{U}}, where U<:Number. δx, dδx and _δv are auxiliary arrays of type Vector{TaylorN{Taylor1{U}}} to avoid allocations. Optionally, the user may provide a Jacobian function jacobianfunc! to compute jac. Otherwise, jac is computed via automatic differentiation using TaylorSeries.jl.

source
TaylorIntegration.stepsize — Method
stepsize(x, epsilon) -> h

Returns a maximum time-step for a the Taylor expansion x using a prescribed absolute tolerance epsilon and the last two Taylor coefficients of (each component of) x.

Note that x is of type Taylor1{U} or Vector{Taylor1{U}}, including also the cases Taylor1{TaylorN{U}} and Vector{Taylor1{TaylorN{U}}}.

Depending of eltype(x), i.e., U<:Number, it may be necessary to overload stepsize, specializing it on the type U, to avoid type instabilities.

source
TaylorIntegration.surfacecrossing — Method
surfacecrossing(g_old, g_now, eventorder::Int)

Detect if the solution crossed a root of event function g. g_old represents the last-before-current value of event function g, and g_now represents the current one; these are Tuple{Bool,Taylor1{T}}s. eventorder is the order of the derivative of the event function g whose root we are trying to find. Returns true if the constant terms of g_old[2] and g_now[2] have different signs (i.e., if one is positive and the other one is negative). Otherwise, if g_old[2] and g_now[2] have the same sign or if the first component of either of them is false, then it returns false.

source
TaylorIntegration.taylorstep! — Method
taylorstep!(f, t, x, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt
-taylorstep!(f!, t, x, dx, xaux, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt

One-step Taylor integration for the one-dependent variable ODE $\dot{x}=dx/dt=f(x, p, t)$ with initial conditions $x(t_0)=x_0$. Returns the time-step δt of the actual integration carried out (δt is positive).

Here, f is the function defining the RHS of the ODE (see taylorinteg), t is the independent variable, x contains the Taylor expansion of the dependent variable, order is the degree used for the Taylor1 polynomials during the integration abstol is the absolute tolerance used to determine the time step of the integration, and params are the parameters entering the ODE functions. For several variables, dx and xaux, both of the same type as x, are needed to save allocations. Finally, parse_eqs is a switch to force not using (parse_eqs=false) the specialized method of jetcoeffs! created with @taylorize; the default is true (parse the equations). Finally, parse_eqs is a switch to force not using (parse_eqs=false) the specialized method of jetcoeffs! created with @taylorize; the default is true (parse the equations).

source

Index

+ newtoniter, nevents) -> nevents

Internal root-finding subroutine, based on Newton-Raphson process. If there is a crossing, then the crossing data is stored in tvS, xvS and gvS and nevents, the number of events/crossings, is updated. Here t is a Taylor1 polynomial which represents the independent variable; x is an array of Taylor1 variables which represent the vector of dependent variables; dx is an array of Taylor1 variables which represent the LHS of the ODE; g_tupl_old is the last-before-current value returned by event function g and g_tupl is the current one; eventorder is the order of the derivative of g whose roots the user is interested in finding; tvS stores the surface-crossing instants; xvS stores the value of the solution at each of the crossings; gvS stores the values of the event function g (or its eventorder-th derivative) at each of the crossings; t0 is the current time; δt_old is the last time-step size; x_dx, x_dx_val, g_dg, g_dg_val are auxiliary variables; nrabstol is the Newton-Raphson process tolerance; newtoniter is the maximum allowed number of Newton-Raphson iteration; nevents is the current number of detected events/crossings.

source
TaylorIntegration.inbookkeeping — Method

inbookkeeping(v, bkkeep::BookKeeping)

Checks if v is declared in bkkeep, considering the d_indx, v_newvars and v_arraydecl fields.

source
TaylorIntegration.jetcoeffs! — Method
jetcoeffs!(eqsdiff!::Function, t, x, dx, xaux, params)

Mutates x in-place using the recursion relation of the derivatives obtained from the differential equations $\dot{x}=dx/dt=f(x, p, t)$.

eqsdiff! is the function defining the RHS of the ODE, x contains the Taylor1 expansion of the dependent variables and t is the independent variable, and params are the parameters appearing on the function defining the differential equation. See taylorinteg for examples and convention for eqsdiff. Note that x is of type Vector{Taylor1{U}} where U<:Number; t is of type Taylor1{T} where T<:Real. In this case, two auxiliary containers dx and xaux (both of the same type as x) are needed to avoid allocations.

Initially, x contains only the 0-th order Taylor coefficient of the current system state (the initial conditions), and jetcoeffs! computes recursively the high-order derivates back into x.

source
TaylorIntegration.jetcoeffs! — Method
jetcoeffs!(eqsdiff::Function, t, x, params)

Returns an updated x using the recursion relation of the derivatives obtained from the differential equations $\dot{x}=dx/dt=f(x, p, t)$.

eqsdiff is the function defining the RHS of the ODE, x contains the Taylor1 expansion of the dependent variable(s) and t is the independent variable, and params are the parameters appearing on the function defining the differential equation. See taylorinteg for examples and convention for eqsdiff. Note that x is of type Taylor1{U} where U<:Number; t is of type Taylor1{T} where T<:Real.

Initially, x contains only the 0-th order Taylor coefficient of the current system state (the initial conditions), and jetcoeffs! computes recursively the high-order derivates back into x.

source
TaylorIntegration.lyap_jetcoeffs! — Method
lyap_jetcoeffs!(t, x, dx, jac, varsaux)

Similar to jetcoeffs! for the calculation of the Lyapunov spectrum. Updates only the elements of x which correspond to the solution of the 1st-order variational equations $\dot{\xi}=J \cdot \xi$, where $J$ is the Jacobian matrix, i.e., the linearization of the equations of motion. jac is the Taylor expansion of $J$ wrt the independent variable, around the current initial condition. varsaux is an auxiliary array of type Array{eltype(jac),3} to avoid allocations. Calling this method assumes that jac has been computed previously using stabilitymatrix!.

source
TaylorIntegration.lyap_taylorstep! — Method
lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, t0, t1, order, abstol, _δv, varsaux, params[, jacobianfunc!])

Similar to taylorstep! for the calculation of the Lyapunov spectrum. jac is the Taylor expansion (wrt the independent variable) of the linearization of the equations of motion, i.e, the Jacobian. xaux, δx, dδx, varsaux and _δv are auxiliary vectors, and params define the parameters of the ODEs. Optionally, the user may provide a Jacobian function jacobianfunc! to compute jac. Otherwise, jac is computed via automatic differentiation using TaylorSeries.jl.

source
TaylorIntegration.nrconvergencecriterion — Method
nrconvergencecriterion(g_val, nrabstol::T, nriter::Int, newtoniter::Int) where {T<:Real}

A rudimentary convergence criterion for the Newton-Raphson root-finding process. g_val may be either a Real, Taylor1{T} or a TaylorN{T}, where T<:Real. Returns true if: 1) the absolute value of g_val, the value of the event function g evaluated at the current estimated root by the Newton-Raphson process, is less than the nrabstol tolerance; and 2) the number of iterations nriter of the Newton-Raphson process is less than the maximum allowed number of iterations, newtoniter; otherwise, returns false.

source
TaylorIntegration.stabilitymatrix! — Method
stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, _δv, params[, jacobianfunc!=nothing])

Updates the matrix jac::Matrix{Taylor1{U}} (linearized equations of motion) computed from the equations of motion (eqsdiff!), at time t at x; x is of type Vector{Taylor1{U}}, where U<:Number. δx, dδx and _δv are auxiliary arrays of type Vector{TaylorN{Taylor1{U}}} to avoid allocations. Optionally, the user may provide a Jacobian function jacobianfunc! to compute jac. Otherwise, jac is computed via automatic differentiation using TaylorSeries.jl.

source
TaylorIntegration.stepsize — Method
stepsize(x, epsilon) -> h

Returns a maximum time-step for a the Taylor expansion x using a prescribed absolute tolerance epsilon and the last two Taylor coefficients of (each component of) x.

Note that x is of type Taylor1{U} or Vector{Taylor1{U}}, including also the cases Taylor1{TaylorN{U}} and Vector{Taylor1{TaylorN{U}}}.

Depending of eltype(x), i.e., U<:Number, it may be necessary to overload stepsize, specializing it on the type U, to avoid type instabilities.

source
TaylorIntegration.surfacecrossing — Method
surfacecrossing(g_old, g_now, eventorder::Int)

Detect if the solution crossed a root of event function g. g_old represents the last-before-current value of event function g, and g_now represents the current one; these are Tuple{Bool,Taylor1{T}}s. eventorder is the order of the derivative of the event function g whose root we are trying to find. Returns true if the constant terms of g_old[2] and g_now[2] have different signs (i.e., if one is positive and the other one is negative). Otherwise, if g_old[2] and g_now[2] have the same sign or if the first component of either of them is false, then it returns false.

source
TaylorIntegration.taylorstep! — Method
taylorstep!(f, t, x, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt
+taylorstep!(f!, t, x, dx, xaux, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt

One-step Taylor integration for the one-dependent variable ODE $\dot{x}=dx/dt=f(x, p, t)$ with initial conditions $x(t_0)=x_0$. Returns the time-step δt of the actual integration carried out (δt is positive).

Here, f is the function defining the RHS of the ODE (see taylorinteg), t is the independent variable, x contains the Taylor expansion of the dependent variable, order is the degree used for the Taylor1 polynomials during the integration abstol is the absolute tolerance used to determine the time step of the integration, and params are the parameters entering the ODE functions. For several variables, dx and xaux, both of the same type as x, are needed to save allocations. Finally, parse_eqs is a switch to force not using (parse_eqs=false) the specialized method of jetcoeffs! created with @taylorize; the default is true (parse the equations). Finally, parse_eqs is a switch to force not using (parse_eqs=false) the specialized method of jetcoeffs! created with @taylorize; the default is true (parse the equations).

source

Index

diff --git a/previews/PR194/common/ddf31eee.svg b/previews/PR194/common/1d00d8fd.svg similarity index 99% rename from previews/PR194/common/ddf31eee.svg rename to previews/PR194/common/1d00d8fd.svg index 147333404..261834e8b 100644 --- a/previews/PR194/common/ddf31eee.svg +++ b/previews/PR194/common/1d00d8fd.svg @@ -1,124 +1,124 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/common/dde1f17e.svg b/previews/PR194/common/4b316cee.svg similarity index 99% rename from previews/PR194/common/dde1f17e.svg rename to previews/PR194/common/4b316cee.svg index 396dc2bea..fed70c8c0 100644 --- a/previews/PR194/common/dde1f17e.svg +++ b/previews/PR194/common/4b316cee.svg @@ -1,45 +1,45 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/common/02d9a208.svg b/previews/PR194/common/6760c6a8.svg similarity index 99% rename from previews/PR194/common/02d9a208.svg rename to previews/PR194/common/6760c6a8.svg index a5fba4208..6b5e0f5a4 100644 --- a/previews/PR194/common/02d9a208.svg +++ b/previews/PR194/common/6760c6a8.svg @@ -1,45 +1,45 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/common/80c901a3.svg b/previews/PR194/common/77e0273c.svg similarity index 99% rename from previews/PR194/common/80c901a3.svg rename to previews/PR194/common/77e0273c.svg index e7b99a385..d3b531c37 100644 --- a/previews/PR194/common/80c901a3.svg +++ b/previews/PR194/common/77e0273c.svg @@ -1,46 +1,46 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/common/4d02b486.svg b/previews/PR194/common/cbc9350b.svg similarity index 92% rename from previews/PR194/common/4d02b486.svg rename to previews/PR194/common/cbc9350b.svg index c49c5408a..fbc1b4423 100644 --- a/previews/PR194/common/4d02b486.svg +++ b/previews/PR194/common/cbc9350b.svg @@ -1,49 +1,49 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/common/index.html b/previews/PR194/common/index.html index 5da553a4e..b6d24715f 100644 --- a/previews/PR194/common/index.html +++ b/previews/PR194/common/index.html @@ -32,7 +32,7 @@ ylims!(-1.7, -1.45) xlabel!("x") ylabel!("J") -title!("Zero-velocity curves (x-axis projection)")Example block output

Notice that the maxima in the plot correspond to the Lagrangian points $L_1$, $L_2$ and $L_3$; below we shall concentrate in the value $J_0 = -1.58$.

J0 = -1.58

We define a function py!, which depends on the initial condition $q_0 = (x_0, 0, 0, p_{y,0})$ and the Jacobi constant value $J_0$, such that it computes an adequate value $p_{y,0}$ for which we have $H(q_0)=J_0$ and updates (in-place) the initial condition accordingly.

function py!(q0, J0)
+title!("Zero-velocity curves (x-axis projection)")
Example block output

Notice that the maxima in the plot correspond to the Lagrangian points $L_1$, $L_2$ and $L_3$; below we shall concentrate in the value $J_0 = -1.58$.

J0 = -1.58

We define a function py!, which depends on the initial condition $q_0 = (x_0, 0, 0, p_{y,0})$ and the Jacobi constant value $J_0$, such that it computes an adequate value $p_{y,0}$ for which we have $H(q_0)=J_0$ and updates (in-place) the initial condition accordingly.

function py!(q0, J0)
     @assert iszero(q0[2]) && iszero(q0[3]) # q0[2] and q0[3] have to be equal to zero
     q0[4] = q0[1] + sqrt( q0[1]^2-2( V(q0[1], q0[2])-J0 ) )
     nothing
@@ -145,39 +145,39 @@
 xlims!(-1+μ-0.2, 1+μ+0.2)
 ylims!(-0.8, 0.8)
 xlabel!("x")
-ylabel!("y")
Example block output

Note that the orbit obtained displays the expected dynamics: the test particle explores the regions surrounding both primaries, located at the red dots, without escaping to infinity. For comparison, we now plot the orbit corresponding to the solution obtained with the Vern9() integration; note that the scales are identical.

plot(solV, vars=(1, 2), linewidth=1, fmt = :png)
+ylabel!("y")
Example block output

Note that the orbit obtained displays the expected dynamics: the test particle explores the regions surrounding both primaries, located at the red dots, without escaping to infinity. For comparison, we now plot the orbit corresponding to the solution obtained with the Vern9() integration; note that the scales are identical.

plot(solV, vars=(1, 2), linewidth=1, fmt = :png)
 scatter!([μ, -1+μ], [0,0], leg=false) # positions of the primaries
 xlims!(-1+μ-0.2, 1+μ+0.2)
 ylims!(-0.8, 0.8)
 xlabel!("x")
-ylabel!("y")
Example block output

We note that both orbits display the same qualitative features, and also some differences. For instance, the TaylorMethod(25) solution gets closer to the primary than that the Vern9(). We can obtain a quantitative comparison of the validity of both integrations through the preservation of the Jacobi constant:

ET = H.(solT.u)
+ylabel!("y")
Example block output

We note that both orbits display the same qualitative features, and also some differences. For instance, the TaylorMethod(25) solution gets closer to the primary than that the Vern9(). We can obtain a quantitative comparison of the validity of both integrations through the preservation of the Jacobi constant:

ET = H.(solT.u)
 EV = H.(solV.u)
 δET = ET .- J0
 δEV = EV .- J0

We plot first the value of the Jacobi constant as function of time.

plot(solT.t, H.(solT.u), label="TaylorMethod(25)", fmt = :png, yformatter = :plain)
 plot!(solV.t, H.(solV.u), label="Vern9()")
 xlabel!("t")
-ylabel!("H")
Example block output

In the scale shown we observe that, while both solutions display a preservation of the Jacobi constant to a certain degree, the Vern9() solution suffers sudden jumps during the integration.

We now plot, in log scale, the abs of the absolute error in the Jacobi constant as a function of time, for both solutions:

plot(solT.t, abs.(δET), yscale=:log10, label="TaylorMethod(25)", legend=:topleft, fmt = :png, yformatter = :plain)
+ylabel!("H")
Example block output

In the scale shown we observe that, while both solutions display a preservation of the Jacobi constant to a certain degree, the Vern9() solution suffers sudden jumps during the integration.

We now plot, in log scale, the abs of the absolute error in the Jacobi constant as a function of time, for both solutions:

plot(solT.t, abs.(δET), yscale=:log10, label="TaylorMethod(25)", legend=:topleft, fmt = :png, yformatter = :plain)
 plot!(solV.t, abs.(δEV), label="Vern9()")
 ylims!(10^-16, 10^-10)
 xlabel!("t")
-ylabel!("dE")
Example block output

We notice that the Jacobi constant absolute error for the TaylorMethod(25) solution remains bounded below $10^{-13}$ throughout the integration. While the Vern9() solution at the end of the integration time has reached a similar value, it displays a larger Jacobi constant error earlier in time.

Finally, we comment on the time spent by each integration.

using BenchmarkTools
+ylabel!("dE")
Example block output

We notice that the Jacobi constant absolute error for the TaylorMethod(25) solution remains bounded below $10^{-13}$ throughout the integration. While the Vern9() solution at the end of the integration time has reached a similar value, it displays a larger Jacobi constant error earlier in time.

Finally, we comment on the time spent by each integration.

using BenchmarkTools
 bT = @benchmark solve($prob, $(TaylorMethod(25)), abstol=1e-15)
-bV = @benchmark solve($prob, $(Vern9()), abstol=1e-15, reltol=1e-15)
bT # TaylorMethod(25) benchmark
BenchmarkTools.Trial: 51 samples with 1 evaluation.
- Range (min … max):  95.620 ms … 112.621 ms  ┊ GC (min … max): 0.00% … 0.00%
- Time  (median):     98.423 ms               ┊ GC (median):    0.00%
- Time  (mean ± σ):   98.357 ms ±   2.482 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%
+bV = @benchmark solve($prob, $(Vern9()), abstol=1e-15, reltol=1e-15)
bT # TaylorMethod(25) benchmark
BenchmarkTools.Trial: 52 samples with 1 evaluation.
+ Range (min … max):  94.465 ms … 102.031 ms  ┊ GC (min … max): 0.00% … 0.00%
+ Time  (median):     97.603 ms               ┊ GC (median):    0.00%
+ Time  (mean ± σ):   97.265 ms ±   1.408 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%
 
-   ▂  ▂                       ▂ █▂█     ▅                       
-  ███▁█▁▅▁▅▁▁▁▅▁▅▁▁▁▁▁▁▁▁▁▁▅▁▁█████▁█▅▅██▅▁▁▁▅▅▁▁▁▁▁▅▁▅▁▅▁▁▁▁▅ ▁
-  95.6 ms         Histogram: frequency by time          101 ms <
+                ▁                 ▁▁▁█ ▁  ▁                     
+  ▄▄▁▄▁▇▇▁▁▁▁▄▄▁█▁▄▁▁▁▁▁▁▁▁▁▁▁▇▁▄▄████▇█▇▄█▇▇▁▄▁▁▄▄▁▁▁▁▁▁▁▁▁▁▄ ▁
+  94.5 ms         Histogram: frequency by time         99.9 ms <
 
  Memory estimate: 7.76 MiB, allocs estimate: 78630.
bV # Vern9 benchmark
BenchmarkTools.Trial: 19 samples with 1 evaluation.
- Range (min … max):  214.017 ms … 651.980 ms  ┊ GC (min … max):  0.00% … 67.10%
- Time  (median):     233.577 ms               ┊ GC (median):     0.00%
- Time  (mean ± σ):   272.096 ms ±  97.655 ms  ┊ GC (mean ± σ):  16.62% ± 18.00%
+ Range (min … max):  215.173 ms … 661.018 ms  ┊ GC (min … max):  0.00% … 67.71%
+ Time  (median):     233.658 ms               ┊ GC (median):     0.00%
+ Time  (mean ± σ):   272.603 ms ±  99.426 ms  ┊ GC (mean ± σ):  16.72% ± 18.06%
 
-    █                                                            
-  ▄▁█▅▁▄▄▁▁▅▁▁▁▁▁▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▁
-  214 ms           Histogram: frequency by time          652 ms <
+    █                                                            
+  ▃▁█▁▃▅▁▁▃▃▁▁▁▁▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
+  215 ms           Histogram: frequency by time          661 ms <
 
- Memory estimate: 130.11 MiB, allocs estimate: 1319694.

We notice in this setup, where the TaylorMethod(25) and the Vern9() integrations perform similarly in terms of accuracy, the former performs better in terms of runtime.

We can tune the abstol and reltol for the Vern9() method we so that performance is similar. Such situation has an accuracy cost, which then makes TaylorIntegration a sensible alternative for high-accuracy integrations of non-stiff ODEs in some cases; see [2].

Finally, as mentioned above, a crucial way in which TaylorIntegration provides high accuracy at competitive speeds is through the use of the @taylorize macro; see this section for details. Currently, TaylorIntegration supports the use of @taylorize via the common interface with DifferentialEquations only for in-place ODEProblem.

References and notes

[1] Murray, Carl D., Stanley F. Dermott. Solar System dynamics. Cambridge University Press, 1999.

[2] SciMLBenchmarks.jl/DynamicalODE

+ Memory estimate: 130.11 MiB, allocs estimate: 1319694.

We notice in this setup, where the TaylorMethod(25) and the Vern9() integrations perform similarly in terms of accuracy, the former performs better in terms of runtime.

We can tune the abstol and reltol for the Vern9() method we so that performance is similar. Such situation has an accuracy cost, which then makes TaylorIntegration a sensible alternative for high-accuracy integrations of non-stiff ODEs in some cases; see [2].

Finally, as mentioned above, a crucial way in which TaylorIntegration provides high accuracy at competitive speeds is through the use of the @taylorize macro; see this section for details. Currently, TaylorIntegration supports the use of @taylorize via the common interface with DifferentialEquations only for in-place ODEProblem.

References and notes

[1] Murray, Carl D., Stanley F. Dermott. Solar System dynamics. Cambridge University Press, 1999.

[2] SciMLBenchmarks.jl/DynamicalODE

diff --git a/previews/PR194/index.html b/previews/PR194/index.html index 2115d18c2..b89240b56 100644 --- a/previews/PR194/index.html +++ b/previews/PR194/index.html @@ -2,4 +2,4 @@ Home · TaylorIntegration.jl
GitHub

TaylorIntegration.jl

ODE integration using Taylor's method, and more, in Julia.

Authors

  • Jorge A. Pérez, Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México (UNAM)
  • Luis Benet, Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México (UNAM)

License

TaylorIntegration is licensed under the MIT "Expat" license.

Installation

TaylorIntegration.jl is a registered package, and is simply installed by running

pkg> add TaylorIntegration

Supporting and citing

This package is developed as part of academic research. If you would like to help supporting it, please star the repository as such metrics may help us secure funding. If you use this software, we would be grateful if you could cite our work as follows (Bibtex entry can be found here):

J.A. Pérez-Hernández and L. Benet
 TaylorIntegration.jl: Taylor Integration in Julia
 https://github.com/PerezHz/TaylorIntegration.jl
-DOI:[10.5281/zenodo.2562352](https://doi.org/10.5281/zenodo.2562352)

Acknowledgments

We acknowledge financial support from the DGAPA-PAPIIT (UNAM) grant IG-100616. LB acknowledges support through a Cátedra Marcos Moshinsky (2013).

+DOI:[10.5281/zenodo.2562352](https://doi.org/10.5281/zenodo.2562352)

Acknowledgments

We acknowledge financial support from the DGAPA-PAPIIT (UNAM) grant IG-100616. LB acknowledges support through a Cátedra Marcos Moshinsky (2013).

diff --git a/previews/PR194/jet_transport/index.html b/previews/PR194/jet_transport/index.html index 129f7831f..01d0c53ca 100644 --- a/previews/PR194/jet_transport/index.html +++ b/previews/PR194/jet_transport/index.html @@ -41,4 +41,4 @@ \delta y \end{array} \right). -\end{aligned}\]

The first terms in the expressions for $\mathbf{x}_1$ and $\mathbf{x}_2$ above correspond to the result of an Euler integration step using the initial conditions only. The other terms are the (linear) corrections which involve the small displacements $\delta x$ and $\delta y$.

In general, for differential equations involving non-linear terms, the resulting expansions in $\delta x$ and $\delta y$ will reflect aspects of the non-linearities of the ODEs. Clearly, jet transport techniques allow to address stability properties beyond the linear case, though memory constraints may play a role. See this example illustrating the implementation for the simple pendulum, and this one illustrating the construction of a Poincaré map with Jet transport techniques.

References

[1] D. Pérez-Palau, Josep J. Masdemont, Gerard Gómez, 2015, Celest. Mech. Dyn. Astron. 123, 239.

+\end{aligned}\]

The first terms in the expressions for $\mathbf{x}_1$ and $\mathbf{x}_2$ above correspond to the result of an Euler integration step using the initial conditions only. The other terms are the (linear) corrections which involve the small displacements $\delta x$ and $\delta y$.

In general, for differential equations involving non-linear terms, the resulting expansions in $\delta x$ and $\delta y$ will reflect aspects of the non-linearities of the ODEs. Clearly, jet transport techniques allow to address stability properties beyond the linear case, though memory constraints may play a role. See this example illustrating the implementation for the simple pendulum, and this one illustrating the construction of a Poincaré map with Jet transport techniques.

References

[1] D. Pérez-Palau, Josep J. Masdemont, Gerard Gómez, 2015, Celest. Mech. Dyn. Astron. 123, 239.

diff --git a/previews/PR194/kepler/2cf8fe75.svg b/previews/PR194/kepler/5381b884.svg similarity index 99% rename from previews/PR194/kepler/2cf8fe75.svg rename to previews/PR194/kepler/5381b884.svg index 382f1a474..2ddab48ff 100644 --- a/previews/PR194/kepler/2cf8fe75.svg +++ b/previews/PR194/kepler/5381b884.svg @@ -1,48 +1,48 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/kepler/c6834815.svg b/previews/PR194/kepler/de4a01a2.svg similarity index 99% rename from previews/PR194/kepler/c6834815.svg rename to previews/PR194/kepler/de4a01a2.svg index f1f4980a0..0b3f021f7 100644 --- a/previews/PR194/kepler/c6834815.svg +++ b/previews/PR194/kepler/de4a01a2.svg @@ -1,50 +1,50 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/kepler/3bcc78c4.svg b/previews/PR194/kepler/faec9260.svg similarity index 97% rename from previews/PR194/kepler/3bcc78c4.svg rename to previews/PR194/kepler/faec9260.svg index 9e4e8356c..d3bc5b417 100644 --- a/previews/PR194/kepler/3bcc78c4.svg +++ b/previews/PR194/kepler/faec9260.svg @@ -1,46 +1,46 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/kepler/index.html b/previews/PR194/kepler/index.html index 55e10d056..e0dab60d6 100644 --- a/previews/PR194/kepler/index.html +++ b/previews/PR194/kepler/index.html @@ -39,7 +39,7 @@ scatter!([0], [0], shape=:circle, ms=5) xaxis!("x", (-2.0, 0.5)) yaxis!("y", (-1.0, 1.0)) -title!("Fig. 1")Example block output

The following functions allow us to calculate the energy and angular momentum using cartesian coordinates.

function energy( x, y, vx, vy )
+title!("Fig. 1")
Example block output

The following functions allow us to calculate the energy and angular momentum using cartesian coordinates.

function energy( x, y, vx, vy )
     kinetic = 0.5 * (vx*vx + vy*vy)
     r = sqrt( x*x + y*y)
     potential = - mu * mass / r
@@ -50,9 +50,9 @@
 plot(t[1:3:end], δE[1:3:end])
 xlabel!("t")
 ylabel!("dE")
-title!("Fig. 2")
Example block output
lz0 = lz(q0...)
+title!("Fig. 2")
Example block output
lz0 = lz(q0...)
 δlz = (lz.(x,y,vx,vy) .- lz0) ./ eps(lz0)
 plot(t[1:3:end], δlz[1:3:end])
 xlabel!("t")
 ylabel!("dlz")
-title!("Fig. 3")
Example block output

These errors are reminiscent of random walks.

The maximum absolute errors of the energy and angular momentum are

maximum( abs.(energy.(x,y,vx,vy) .- e0) ), maximum( abs.(lz.(x,y,vx,vy) .- lz0) )
(9.947598300641403e-14, 2.808864252301646e-14)
+title!("Fig. 3")Example block output

These errors are reminiscent of random walks.

The maximum absolute errors of the energy and angular momentum are

maximum( abs.(energy.(x,y,vx,vy) .- e0) ), maximum( abs.(lz.(x,y,vx,vy) .- lz0) )
(9.947598300641403e-14, 2.808864252301646e-14)
diff --git a/previews/PR194/lorenz_lyapunov/9cda1cfe.svg b/previews/PR194/lorenz_lyapunov/80fb39eb.svg similarity index 97% rename from previews/PR194/lorenz_lyapunov/9cda1cfe.svg rename to previews/PR194/lorenz_lyapunov/80fb39eb.svg index 2b3a24cea..fc70918cd 100644 --- a/previews/PR194/lorenz_lyapunov/9cda1cfe.svg +++ b/previews/PR194/lorenz_lyapunov/80fb39eb.svg @@ -1,50 +1,50 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/lorenz_lyapunov/368faded.svg b/previews/PR194/lorenz_lyapunov/ef7874ee.svg similarity index 94% rename from previews/PR194/lorenz_lyapunov/368faded.svg rename to previews/PR194/lorenz_lyapunov/ef7874ee.svg index 38e17274c..e5d452493 100644 --- a/previews/PR194/lorenz_lyapunov/368faded.svg +++ b/previews/PR194/lorenz_lyapunov/ef7874ee.svg @@ -1,62 +1,62 @@ - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/lorenz_lyapunov/index.html b/previews/PR194/lorenz_lyapunov/index.html index 1b3e384ab..e5e402ed3 100644 --- a/previews/PR194/lorenz_lyapunov/index.html +++ b/previews/PR194/lorenz_lyapunov/index.html @@ -39,10 +39,10 @@ nothing end
Note

We use of zero(x[1]) in the function lorenz_jac! when the RHS consists of a numeric value; this is needed to allow the proper promotion of the variables to carry out Taylor's method.

We can actually check the consistency of lorenz_jac! with the computation of the jacobian using automatic differentiation techniques. Below we use the initial conditions x0, but it is easy to generalize this.

lorenz_jac!(jjac, x0, params, t0)  # update the matrix `jjac` using Jacobian provided by the user
 TaylorSeries.jacobian(dx0TN) == jjac    # `dx0TN` is obtained via automatic differentiation
true

Now, we are ready to perform the integration using lyap_taylorinteg function, which integrates the 1st variational equations and uses Oseledets' theorem. The expansion order will be $28$ and the local absolute tolerance will be $10^{-20}$. lyap_taylorinteg will return three arrays: one with the evaluation times, one with the values of the dependent variables (at the time of evaluation), and another one with the values of the Lyapunov spectrum.

We first carry out the integration computing internally the Jacobian

tv, xv, λv = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params; maxsteps=2000000);

Now, the integration is obtained exploiting lorenz_jac!:

tv_, xv_, λv_ = lyap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20, params, lorenz_jac!; maxsteps=2000000);

In terms of performance the second method is about ~50% faster than the first.

We check the consistency of the orbits computed by the two methods:

tv == tv_, xv == xv_, λv == λv_
(true, true, true)

As mentioned above, a more subtle check is related to the fact that the trace of the Jacobian is constant in time, which must coincide with the sum of all Lyapunov exponents. Using its initial value lorenztr, we compare it with the final Lyapunov exponents of the computation, and obtain

sum(λv[end,:]) ≈ lorenztr, sum(λv_[end,:]) ≈ lorenztr, sum(λv[end,:]) ≈ sum(λv_[end,:])
(true, true, true)

Above we checked the approximate equality; we now show that the relative error is quite small and comparable with the local machine epsilon value around lorenztr:

abs(sum(λv[end,:])/lorenztr - 1), abs(sum(λv_[end,:])/lorenztr - 1), eps(lorenztr)
(1.1102230246251565e-15, 1.1102230246251565e-15, 3.552713678800501e-15)

Therefore, the numerical error is dominated by roundoff errors in the floating point arithmetic of the integration. We will now proceed to plot our results. First, we plot Lorenz attractor in phase space

using Plots
-plot(xv[:,1], xv[:,2], xv[:,3], leg=false)
Example block output

We display now the Lyapunov exponents as a function of time:

using Plots
+plot(xv[:,1], xv[:,2], xv[:,3], leg=false)
Example block output

We display now the Lyapunov exponents as a function of time:

using Plots
 plot(tv, λv[:,1], label="L_1", legend=:right)
 plot!(tv, λv[:,2], label="L_2")
 plot!(tv, λv[:,3], label="L_3")
 xlabel!("time")
 ylabel!("L_i, i=1,2,3")
-title!("Lyapunov exponents vs time")
Example block output

This plot shows that the calculation of the Lyapunov exponents has converged.

+title!("Lyapunov exponents vs time")Example block output

This plot shows that the calculation of the Lyapunov exponents has converged.

diff --git a/previews/PR194/lyapunov_spectrum/index.html b/previews/PR194/lyapunov_spectrum/index.html index 60e3fda4c..d720c9d3f 100644 --- a/previews/PR194/lyapunov_spectrum/index.html +++ b/previews/PR194/lyapunov_spectrum/index.html @@ -1,2 +1,2 @@ -Lyapunov spectrum · TaylorIntegration.jl
GitHub

Lyapunov spectrum

Here we describe the background of the Lyapunov spectra computations in TaylorIntegration.jl. Our implementation follows the numerical method of Benettin et al. [1], [2], which itself is based on Oseledet's multiplicative ergodic theorem [3]. Namely, simultaneously to the integration of the equations of motion, we integrate the 1st-order variational equations associated to them.

In general, given a dynamical system defined by the equations of motion

\[\dot{x} = f(t, x),\tag{1}\]

along with the initial condition $x(t_0) = x_0$, then the first-order variational equations associated to this system are

\[\dot{\xi} = (\operatorname{D}f)(x(t))\cdot \xi,\tag{2}\]

where $(\operatorname{D}f)(x(t))$ is the Jacobian of the function $f$ with respect to the dependent variable $x$, evaluated at time $t$, for a given solution $x(t)$ to the equations of motion. The variable $\xi$ denotes a matrix, whose initial condition is $\xi(t_0) = \mathbb{1}_{n}$, the $n\times n$ identity matrix, where $n$ is the degrees of freedom or number of dependent variables $x$.

For the actual computation of the Lyapunov spectrum, we proceed as follows. During the simultaneous numerical integration of the equations of motion and the variational equations, at fixed time intervals $t_k = k\cdot \Delta t$, $k = 1, 2, \ldots$ we perform a $QR$ decomposition over $\xi(t_k)$, the solution of the variational equations at time $t_k$. That is, we factorize $\xi(t_k)$ as $\xi(t_k)=Q_k\cdot R_k$, where $Q_k$ is an orthogonal $n\times n$ matrix and $R_k$ is an upper triangular $n\times n$ matrix with positive diagonal elements. The diagonal elements $R_{ii,k}$ are the growth factors from which the $l$-th Lyapunov exponent is computed at time $t_k$

\[\lambda_l = \sum_{m=1}^k \frac{\log (R_{ll,m})}{k\cdot \Delta t}.\tag{3}\]

In turn, the matrix $Q$ is substituted into $\xi(t_k)$ as the new (scaled) initial condition.

The equations of motion together with the variational equations are integrated up to time $t_{k+1}$ using Taylor's method. We note that each time step of the integration is determined using the normalized derivatives of $x$ and the tolerance $\epsilon_\textrm{tol}$. This process is repeated until a prescribed $t_\textrm{max}$ is reached.

This example illustrates the computation of the Lyapunov spectrum for the Lorenz system.

References

[1] Benettin G., Galgani L., Giorgilli A., Strelcyn J.M., 1980, Meccanica, 15, 9

[2] Benettin G., Galgani L., Giorgilli A., Strelcyn J.M., 1980, Meccanica, 15, 21

[3] Oseledets V. I., 1968, Trudy Moskovskogo Matematicheskogo Obshchestva, 19, 179

+Lyapunov spectrum · TaylorIntegration.jl
GitHub

Lyapunov spectrum

Here we describe the background of the Lyapunov spectra computations in TaylorIntegration.jl. Our implementation follows the numerical method of Benettin et al. [1], [2], which itself is based on Oseledet's multiplicative ergodic theorem [3]. Namely, simultaneously to the integration of the equations of motion, we integrate the 1st-order variational equations associated to them.

In general, given a dynamical system defined by the equations of motion

\[\dot{x} = f(t, x),\tag{1}\]

along with the initial condition $x(t_0) = x_0$, then the first-order variational equations associated to this system are

\[\dot{\xi} = (\operatorname{D}f)(x(t))\cdot \xi,\tag{2}\]

where $(\operatorname{D}f)(x(t))$ is the Jacobian of the function $f$ with respect to the dependent variable $x$, evaluated at time $t$, for a given solution $x(t)$ to the equations of motion. The variable $\xi$ denotes a matrix, whose initial condition is $\xi(t_0) = \mathbb{1}_{n}$, the $n\times n$ identity matrix, where $n$ is the degrees of freedom or number of dependent variables $x$.

For the actual computation of the Lyapunov spectrum, we proceed as follows. During the simultaneous numerical integration of the equations of motion and the variational equations, at fixed time intervals $t_k = k\cdot \Delta t$, $k = 1, 2, \ldots$ we perform a $QR$ decomposition over $\xi(t_k)$, the solution of the variational equations at time $t_k$. That is, we factorize $\xi(t_k)$ as $\xi(t_k)=Q_k\cdot R_k$, where $Q_k$ is an orthogonal $n\times n$ matrix and $R_k$ is an upper triangular $n\times n$ matrix with positive diagonal elements. The diagonal elements $R_{ii,k}$ are the growth factors from which the $l$-th Lyapunov exponent is computed at time $t_k$

\[\lambda_l = \sum_{m=1}^k \frac{\log (R_{ll,m})}{k\cdot \Delta t}.\tag{3}\]

In turn, the matrix $Q$ is substituted into $\xi(t_k)$ as the new (scaled) initial condition.

The equations of motion together with the variational equations are integrated up to time $t_{k+1}$ using Taylor's method. We note that each time step of the integration is determined using the normalized derivatives of $x$ and the tolerance $\epsilon_\textrm{tol}$. This process is repeated until a prescribed $t_\textrm{max}$ is reached.

This example illustrates the computation of the Lyapunov spectrum for the Lorenz system.

References

[1] Benettin G., Galgani L., Giorgilli A., Strelcyn J.M., 1980, Meccanica, 15, 9

[2] Benettin G., Galgani L., Giorgilli A., Strelcyn J.M., 1980, Meccanica, 15, 21

[3] Oseledets V. I., 1968, Trudy Moskovskogo Matematicheskogo Obshchestva, 19, 179

diff --git a/previews/PR194/pendulum/868f8045.svg b/previews/PR194/pendulum/720c6c0e.svg similarity index 88% rename from previews/PR194/pendulum/868f8045.svg rename to previews/PR194/pendulum/720c6c0e.svg index 1be6649a3..b83f18b08 100644 --- a/previews/PR194/pendulum/868f8045.svg +++ b/previews/PR194/pendulum/720c6c0e.svg @@ -1,140 +1,140 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/pendulum/index.html b/previews/PR194/pendulum/index.html index d77401dd0..d57c8d974 100644 --- a/previews/PR194/pendulum/index.html +++ b/previews/PR194/pendulum/index.html @@ -34,4 +34,4 @@ scatter!( x_nom, p_nom, color=:black, m=(1,2.8,stroke(0)) -)Example block output +)Example block output diff --git a/previews/PR194/poincareanim1.gif b/previews/PR194/poincareanim1.gif index dbbab94be..2c1db416e 100644 Binary files a/previews/PR194/poincareanim1.gif and b/previews/PR194/poincareanim1.gif differ diff --git a/previews/PR194/poincareanim3.gif b/previews/PR194/poincareanim3.gif index b1c8e1e00..06f3e32b4 100644 Binary files a/previews/PR194/poincareanim3.gif and b/previews/PR194/poincareanim3.gif differ diff --git a/previews/PR194/root_finding/index.html b/previews/PR194/root_finding/index.html index 99906ff1c..5add5ab47 100644 --- a/previews/PR194/root_finding/index.html +++ b/previews/PR194/root_finding/index.html @@ -66,7 +66,7 @@ xlabel!("y") ylabel!("py") title!("Hénon-Heiles Poincaré map (21 iterates)") -end
Plots.Animation("/tmp/jl_3GaLsZ", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000012.png", "000013.png", "000014.png", "000015.png", "000016.png", "000017.png", "000018.png", "000019.png", "000020.png", "000021.png"])
gif(poincare_anim1, "poincareanim1.gif", fps = 2)
[ Info: Saved animation to /home/runner/work/TaylorIntegration.jl/TaylorIntegration.jl/docs/build/poincareanim1.gif

Poincaré map for the Hénon Heiles system

Jet transport

Now, we illustrate the use of jet transport techniques in the same example, that is, we propagate a neighborhood around x0, which will be plotted in the Poincaré map. We first define the vector of small increments of the phase space variables, xTN; we fix the maximum order of the polynomial expansion in these variables to be 4. Then, x0TN is the neighborhood in the 4-dimensional phase space around $x0$.

xTN = set_variables("δx δy δpx δpy", numvars=length(x0), order=4)
+end
Plots.Animation("/tmp/jl_uA3TnX", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000012.png", "000013.png", "000014.png", "000015.png", "000016.png", "000017.png", "000018.png", "000019.png", "000020.png", "000021.png"])
gif(poincare_anim1, "poincareanim1.gif", fps = 2)
[ Info: Saved animation to /home/runner/work/TaylorIntegration.jl/TaylorIntegration.jl/docs/build/poincareanim1.gif

Poincaré map for the Hénon Heiles system

Jet transport

Now, we illustrate the use of jet transport techniques in the same example, that is, we propagate a neighborhood around x0, which will be plotted in the Poincaré map. We first define the vector of small increments of the phase space variables, xTN; we fix the maximum order of the polynomial expansion in these variables to be 4. Then, x0TN is the neighborhood in the 4-dimensional phase space around $x0$.

xTN = set_variables("δx δy δpx δpy", numvars=length(x0), order=4)
 x0TN = x0 .+ xTN

As it was shown above, $x0$ belongs to the energy surface $H(x0) = E_0 = 0.1025$; yet, as it was defined above, the set of phase space points denoted by x0TN includes points that belong to other energy surfaces. This can be noticed by computing H(x0TN)

H(x0TN)
 0.1025 + 0.2478237775 δy + 0.24817464176177093 δpx + 0.9533499999999999 δx² + 0.046650000000000025 δy² + 0.5 δpx² + 0.5 δpy² + 1.0 δx² δy - 0.3333333333333333 δy³ + 𝒪(‖x‖⁵)

Clearly, the expression above may contain points whose energy is different from E0. As it was done above, we shall fix the px component of x0TN so all points of the neighborhood are in the same energy surface.

px!(x0TN, E0) # Impose that all variations are on the proper energy shell!
 H(x0TN)
 0.1025 - 5.551115123125783e-17 δy² - 5.551115123125783e-17 δy³ - 1.7763568394002505e-15 δx² δy² - 4.440892098500626e-16 δy⁴ + 8.881784197001252e-16 δy² δpy² + 𝒪(‖x‖⁵)

We notice that the coefficients of all monomials whose order is not zero are very small, and the constant_term is E0.

In order to properly handle this case, we need to extend the definition of g to be useful for Taylor1{TaylorN{T}} vectors.

#specialized method of g for Taylor1{TaylorN{T}}'s
 function g(dx::Array{Taylor1{TaylorN{T}},1}, x::Array{Taylor1{TaylorN{T}},1},
@@ -116,11 +116,11 @@
     xlabel!("y")
     ylabel!("py")
     title!("Poincaré map: 4th-order jet transport vs Monte Carlo")
-end
Plots.Animation("/tmp/jl_x1Wje0", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000012.png", "000013.png", "000014.png", "000015.png", "000016.png", "000017.png", "000018.png", "000019.png", "000020.png", "000021.png"])
gif(poincare_anim2, "poincareanim2.gif", fps = 2)
[ Info: Saved animation to /home/runner/work/TaylorIntegration.jl/TaylorIntegration.jl/docs/build/poincareanim2.gif

Poincaré map: Jet transport vs Monte Carlo

The next animation is the same as before, adapting the scale.

poincare_anim3 = @animate for i=1:21
+end
Plots.Animation("/tmp/jl_Ivb9BU", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000012.png", "000013.png", "000014.png", "000015.png", "000016.png", "000017.png", "000018.png", "000019.png", "000020.png", "000021.png"])
gif(poincare_anim2, "poincareanim2.gif", fps = 2)
[ Info: Saved animation to /home/runner/work/TaylorIntegration.jl/TaylorIntegration.jl/docs/build/poincareanim2.gif

Poincaré map: Jet transport vs Monte Carlo

The next animation is the same as before, adapting the scale.

poincare_anim3 = @animate for i=1:21
     scatter(map(x->x[i,2], xvSv), map(x->x[i,4], xvSv), marker=(:circle, stroke(0)),
         markersize=0.1, label="Monte Carlo")
     plot!(yS[i,:], pS[i,:], width=0.5, label="Jet transport")
     xlabel!("y")
     ylabel!("py")
     title!("Poincaré map: 4th-order jet transport vs Monte Carlo")
-end
Plots.Animation("/tmp/jl_t2IjZ0", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000012.png", "000013.png", "000014.png", "000015.png", "000016.png", "000017.png", "000018.png", "000019.png", "000020.png", "000021.png"])
gif(poincare_anim3, "poincareanim3.gif", fps = 2)
[ Info: Saved animation to /home/runner/work/TaylorIntegration.jl/TaylorIntegration.jl/docs/build/poincareanim3.gif

Poincaré map: Jet transport vs Monte Carlo

+end
Plots.Animation("/tmp/jl_PfrmSh", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000012.png", "000013.png", "000014.png", "000015.png", "000016.png", "000017.png", "000018.png", "000019.png", "000020.png", "000021.png"])
gif(poincare_anim3, "poincareanim3.gif", fps = 2)
[ Info: Saved animation to /home/runner/work/TaylorIntegration.jl/TaylorIntegration.jl/docs/build/poincareanim3.gif

Poincaré map: Jet transport vs Monte Carlo

diff --git a/previews/PR194/simple_example/331b6c7a.svg b/previews/PR194/simple_example/0d30ad0b.svg similarity index 74% rename from previews/PR194/simple_example/331b6c7a.svg rename to previews/PR194/simple_example/0d30ad0b.svg index 0a3789f58..7f86dd76e 100644 --- a/previews/PR194/simple_example/331b6c7a.svg +++ b/previews/PR194/simple_example/0d30ad0b.svg @@ -1,188 +1,188 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/simple_example/d7a6ff25.svg b/previews/PR194/simple_example/ee2cfd3a.svg similarity index 73% rename from previews/PR194/simple_example/d7a6ff25.svg rename to previews/PR194/simple_example/ee2cfd3a.svg index 13ff0e056..40a7d51bc 100644 --- a/previews/PR194/simple_example/d7a6ff25.svg +++ b/previews/PR194/simple_example/ee2cfd3a.svg @@ -1,200 +1,200 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/previews/PR194/simple_example/index.html b/previews/PR194/simple_example/index.html index 0a5ca8809..ba1f0dcd2 100644 --- a/previews/PR194/simple_example/index.html +++ b/previews/PR194/simple_example/index.html @@ -6,12 +6,12 @@ xlabel!("t") ylabel!("log10(x(t))") xlims!(0,0.34) -title!("Fig. 1")Example block output

Clearly, the solution diverges without bound when $t\to t_\textrm{max} = 1/3$, i.e., $x(t)$ approaches infinity in finite time.

Figure 2 shows the relative difference between the numerical and the analytical solution in terms of time.

exactsol(t, x0) = x0 / (1 - x0 * t)
+title!("Fig. 1")
Example block output

Clearly, the solution diverges without bound when $t\to t_\textrm{max} = 1/3$, i.e., $x(t)$ approaches infinity in finite time.

Figure 2 shows the relative difference between the numerical and the analytical solution in terms of time.

exactsol(t, x0) = x0 / (1 - x0 * t)
 δxT = abs.(xT .- exactsol.(tT, 3.0)) ./ exactsol.(tT, 3.0);
 plot(tT[6:end], log10.(δxT[6:end]), shape=:circle)
 xlabel!("t")
 ylabel!("log10(dx(t))")
 xlims!(0, 0.4)
-title!("Fig. 2")
Example block output

To put in perspective how good is the constructed solution, we impose (arbitrarily) a relative accuracy of $10^{-13}$; the time until such accuracy is satisfied is given by:

indx = findfirst(δxT .> 1.0e-13);
+title!("Fig. 2")
Example block output

To put in perspective how good is the constructed solution, we impose (arbitrarily) a relative accuracy of $10^{-13}$; the time until such accuracy is satisfied is given by:

indx = findfirst(δxT .> 1.0e-13);
 esol = exactsol(tT[indx-1],3.0);
-tT[indx-1], esol, eps(esol)
(0.3323590833381362, 1026.4305926915836, 2.2737367544323206e-13)

Note that, the accuracy imposed in terms of the actual value of the exact solution means that the difference of the computed and the exact solution is essentially due to the eps of the computed value.

+tT[indx-1], esol, eps(esol)
(0.3323590833381362, 1026.4305926915836, 2.2737367544323206e-13)

Note that, the accuracy imposed in terms of the actual value of the exact solution means that the difference of the computed and the exact solution is essentially due to the eps of the computed value.

diff --git a/previews/PR194/taylor_method/index.html b/previews/PR194/taylor_method/index.html index 56f145b54..e28701658 100644 --- a/previews/PR194/taylor_method/index.html +++ b/previews/PR194/taylor_method/index.html @@ -1,4 +1,4 @@ Taylor's method · TaylorIntegration.jl
GitHub

Taylor's method

Taylor's integration method is a quite powerful method to integrate ODEs which are smooth enough, allowing to reach a precision comparable to round-off errors per time-step. A high-order Taylor approximation of the solution (dependent variable) is constructed such that the error is quite small. A time-step is constructed which guarantees the validity of the series; this is used to sum up the Taylor expansion to obtain an approximation of the solution at a later time.

The recurrence relation

Let us consider the following

\[\dot{x} = f(t, x),\tag{1}\]

and define the initial value problem with the initial condition $x(t_0) = x(0)$.

We write the solution of this equation as

\[x = x_{[0]} + x_{[1]} (t-t_0) + x_{[2]} (t-t_0)^2 + \cdots + x_{[k]} (t-t_0)^k + \cdots,\tag{2}\]

where the initial condition imposes that $x_{[0]} = x(0)$. Below, we show how to obtain the coefficients $x_{[k]}$ of the Taylor expansion of the solution.

We assume that the Taylor expansion around $t_0$ of $f(t, x(t))$ is known, which we write as

\[f(t, x(t)) = f_{[0]} + f_{[1]} (t-t_0) + f_{[2]} (t-t_0)^2 + \cdots -+ f_{[k]} (t-t_0)^k + \cdots.\tag{3}\]

Here, $f_{[0]}=f(t_0,x_0)$, and the Taylor coefficients $f_{[k]} = f_{[k]}(t_0)$ are the $k$-th normalized derivatives at $t_0$ given by

\[f_{[k]} = \frac{1}{k!} \frac{{\rm d}^k f} {{\rm d} t^k}(t_0).\tag{4}\]

Then, we are assuming that we know how to obtain $f_{[k]}$; these coefficients are obtained using TaylorSeries.jl.

Substituting Eq. (2) in (1), and equating powers of $t-t_0$, we obtain

\[x_{[k+1]} = \frac{f_{[k]}(t_0)}{k+1}, \quad k=0,1,\dots.\tag{5}\]

Therefore, the coefficients of the Taylor expansion (2) are obtained recursively using Eq. (5).

Time step

In the computer, the expansion (2) has to be computed to a finite order. We shall denote by $K$ the order of the series. In principle, the larger the order $K$, the more accurate the obtained solution is.

The theorem of existence and uniqueness of the solution of Eq.~(1) ensures that the Taylor expansion converges. Then, assuming that $K$ is large enough to be within the convergent tail, we introduce the parameter $\epsilon_\textrm{tol} > 0$ to control how large is the last term. The idea is to set this parameter to a small value, usually smaller than the machine-epsilon. Denoting by $h = t_1-t_0$ the time step, then imposing $| x_{[K]} | h^K \le \epsilon_\textrm{tol}$ we obtain

\[h \le \Big(\frac{\epsilon_\textrm{tol}}{| x_{[K]} |}\Big)^{1/K}.\tag{6}\]

Equation (6) represents the maximum time-step which is consistent with $\epsilon_\textrm{tol}$, $K$ and the assumption of being within the convergence tail. Notice that the arguments exposed above simply ensure that $h$ is a maximum time-step, but any other smaller than $h$ can be used since the series is convergent in the open interval $t\in(t_0-h,t_0+h)$.

Finally, from Eq. (2) with (6) we obtain $x(t_1) = x(t_0+h)$, which is again an initial value problem.

++ f_{[k]} (t-t_0)^k + \cdots.\tag{3}\]

Here, $f_{[0]}=f(t_0,x_0)$, and the Taylor coefficients $f_{[k]} = f_{[k]}(t_0)$ are the $k$-th normalized derivatives at $t_0$ given by

\[f_{[k]} = \frac{1}{k!} \frac{{\rm d}^k f} {{\rm d} t^k}(t_0).\tag{4}\]

Then, we are assuming that we know how to obtain $f_{[k]}$; these coefficients are obtained using TaylorSeries.jl.

Substituting Eq. (2) in (1), and equating powers of $t-t_0$, we obtain

\[x_{[k+1]} = \frac{f_{[k]}(t_0)}{k+1}, \quad k=0,1,\dots.\tag{5}\]

Therefore, the coefficients of the Taylor expansion (2) are obtained recursively using Eq. (5).

Time step

In the computer, the expansion (2) has to be computed to a finite order. We shall denote by $K$ the order of the series. In principle, the larger the order $K$, the more accurate the obtained solution is.

The theorem of existence and uniqueness of the solution of Eq.~(1) ensures that the Taylor expansion converges. Then, assuming that $K$ is large enough to be within the convergent tail, we introduce the parameter $\epsilon_\textrm{tol} > 0$ to control how large is the last term. The idea is to set this parameter to a small value, usually smaller than the machine-epsilon. Denoting by $h = t_1-t_0$ the time step, then imposing $| x_{[K]} | h^K \le \epsilon_\textrm{tol}$ we obtain

\[h \le \Big(\frac{\epsilon_\textrm{tol}}{| x_{[K]} |}\Big)^{1/K}.\tag{6}\]

Equation (6) represents the maximum time-step which is consistent with $\epsilon_\textrm{tol}$, $K$ and the assumption of being within the convergence tail. Notice that the arguments exposed above simply ensure that $h$ is a maximum time-step, but any other smaller than $h$ can be used since the series is convergent in the open interval $t\in(t_0-h,t_0+h)$.

Finally, from Eq. (2) with (6) we obtain $x(t_1) = x(t_0+h)$, which is again an initial value problem.

diff --git a/previews/PR194/taylorize/index.html b/previews/PR194/taylorize/index.html index 6e69390f3..6a5846247 100644 --- a/previews/PR194/taylorize/index.html +++ b/previews/PR194/taylorize/index.html @@ -16,7 +16,7 @@ t1, x1 = taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run e1 = @elapsed taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000); all1 = @allocated taylorinteg(pendulumNP!, q0, t0, tf, 25, 1e-20, maxsteps=50000); -e1, all1
(0.89130148, 761690960)

The initial number of methods defined for TaylorIntegration.jetcoeffs! is 2; yet, since @taylorize was used in an example previously, the current number of methods is 3, as explained below.

println(methods(TaylorIntegration.jetcoeffs!)) # default methods
# 3 methods for generic function "jetcoeffs!" from TaylorIntegration:
+e1, all1
(0.895614484, 761690960)

The initial number of methods defined for TaylorIntegration.jetcoeffs! is 2; yet, since @taylorize was used in an example previously, the current number of methods is 3, as explained below.

println(methods(TaylorIntegration.jetcoeffs!)) # default methods
# 3 methods for generic function "jetcoeffs!" from TaylorIntegration:
  [1] jetcoeffs!(::Val{Main.__atexample__named__common.pcr3bp!}, t::Taylor1{_T}, q::AbstractArray{Taylor1{_S}, _N}, dq::AbstractArray{Taylor1{_S}, _N}, param, __ralloc::TaylorIntegration.RetAlloc{Taylor1{_S}}) where {_T<:Real, _S<:Number, _N}
      @ Main.__atexample__named__common none:0
  [2] jetcoeffs!(eqsdiff!::Function, t::Taylor1{T}, x::AbstractArray{Taylor1{U}, N}, dx::AbstractArray{Taylor1{U}, N}, xaux::AbstractArray{Taylor1{U}, N}, params) where {T<:Real, U<:Number, N}
@@ -47,17 +47,17 @@
      @ ~/work/TaylorIntegration.jl/TaylorIntegration.jl/src/integrator.jl:23

We see that there is only one method of pendulum!, and there is a new method (four in total) of TaylorIntegration.jetcoeffs!, whose signature appears in this documentation as Val{Main.pendulum!}. It is a specialized version for the function pendulum! (with some extra information about the module where the function was created). This method is selected internally if it exists (default), exploiting dispatch, when calling taylorinteg or lyap_taylorinteg. In order to integrate using the hard-coded standard (default) method of TaylorIntegration.jetcoeffs! of the integration above, the keyword argument parse_eqs has to be set to false. Similarly, one can check that there exists a new method of TaylorIntegration._allocate_jetcoeffs!.

Now we carry out the integration using the specialized method; note that we use the same instruction as above; the default value for the keyword argument parse_eqs is true, so we may omit it.

t2, x2 = taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000); # warm-up run
 e2 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
 all2 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000);
-e2, all2
(0.055402838, 1210896)

We note the difference in the performance and allocations:

e1/e2, all1/all2
(16.08765023914479, 629.0308663997569)

We can check that both integrations yield the same results.

t1 == t2 && x1 == x2
true

As stated above, in order to allow opting out of using the specialized method created by @taylorize, taylorinteg and lyap_taylorinteg recognize the keyword argument parse_eqs; setting it to false causes the standard method to be used.

taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000, parse_eqs=false); # warm-up run
+e2, all2
(0.05502757, 1210896)

We note the difference in the performance and allocations:

e1/e2, all1/all2
(16.275741123949324, 629.0308663997569)

We can check that both integrations yield the same results.

t1 == t2 && x1 == x2
true

As stated above, in order to allow opting out of using the specialized method created by @taylorize, taylorinteg and lyap_taylorinteg recognize the keyword argument parse_eqs; setting it to false causes the standard method to be used.

taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000, parse_eqs=false); # warm-up run
 
 e3 = @elapsed taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000, parse_eqs=false);
 all3 = @allocated taylorinteg(pendulum!, q0, t0, tf, 25, 1e-20, maxsteps=50000, parse_eqs=false);
-e1/e3, all1/all3
(0.9904659063654225, 1.0000025837317077)

We now illustrate the possibility of exploiting the macro when using TaylorIntegration.jl from DifferentialEquations.jl.

using OrdinaryDiffEq
+e1/e3, all1/all3
(0.9976122354622782, 1.0000025837317077)

We now illustrate the possibility of exploiting the macro when using TaylorIntegration.jl from DifferentialEquations.jl.

using OrdinaryDiffEq
 
 prob = ODEProblem(pendulum!, q0, (t0, tf), nothing) # no parameters
 solT = solve(prob, TaylorMethod(25), abstol=1e-20, parse_eqs=true); # warm-up run
 e4 = @elapsed solve(prob, TaylorMethod(25), abstol=1e-20, parse_eqs=true);
 
-e1/e4
13.854473307908483

Note that there is an additional cost to using solve in comparison with taylorinteg, but still @taylorize yields improved running times.

The speed-up obtained comes from the design of the new (specialized) method of TaylorIntegration.jetcoeffs! as described above: it avoids some allocations and some repeated computations. This is achieved by knowing the specific AST of the function of the ODEs integrated, which is walked through and translated into the actual implementation, where some required auxiliary arrays are created and reused, and the low-level functions defined in TaylorSeries.jl are used. For this, we heavily rely on Espresso.jl and some metaprogramming; we thank Andrei Zhabinski for his help and comments.

The new TaylorIntegration.jetcoeffs! and TaylorIntegration._allocate_jetcoeffs! methods can be inspected by constructing the expression corresponding to the function, and using TaylorIntegration._make_parsed_jetcoeffs:

ex = :(function pendulum!(dx, x, p, t)
+e1/e4
13.94104803449267

Note that there is an additional cost to using solve in comparison with taylorinteg, but still @taylorize yields improved running times.

The speed-up obtained comes from the design of the new (specialized) method of TaylorIntegration.jetcoeffs! as described above: it avoids some allocations and some repeated computations. This is achieved by knowing the specific AST of the function of the ODEs integrated, which is walked through and translated into the actual implementation, where some required auxiliary arrays are created and reused, and the low-level functions defined in TaylorSeries.jl are used. For this, we heavily rely on Espresso.jl and some metaprogramming; we thank Andrei Zhabinski for his help and comments.

The new TaylorIntegration.jetcoeffs! and TaylorIntegration._allocate_jetcoeffs! methods can be inspected by constructing the expression corresponding to the function, and using TaylorIntegration._make_parsed_jetcoeffs:

ex = :(function pendulum!(dx, x, p, t)
     dx[1] = x[2]
     dx[2] = -sin(x[1])
     return dx
@@ -95,4 +95,4 @@
       tmp378 = Taylor1(cos(constant_term(x[1])), order)
       dx[2] = Taylor1(-(constant_term(tmp376)), order)
       return TaylorIntegration.RetAlloc{Taylor1{_S}}([tmp376, tmp378], [Array{Taylor1{_S}, 1}(undef, 0)], [Array{Taylor1{_S}, 2}(undef, 0, 0)], [Array{Taylor1{_S}, 3}(undef, 0, 0, 0)], [Array{Taylor1{_S}, 4}(undef, 0, 0, 0, 0)])
-  end))

The first function has a similar structure as the hard-coded method of TaylorIntegration.jetcoeffs!, but uses low-level functions in TaylorSeries (e.g., sincos! above). Temporary Taylor1 objects as well as declared arrays are allocated once by TaylorIntegration._allocate_jetcoeffs!. More complex functions quickly become very difficult to read. Note that, if necessary, one can further optimize new_ex manually.

Limitations and some advice

The construction of the internal function obtained by using @taylorize is somewhat complicated and limited. Here we list some limitations and provide some advice.

  • Expressions must involve two arguments at most, which requires using parentheses appropriately: For example, res = a+b+c should be written as res = (a+b)+c. This may lead to more parentheses used in compound expressions than would be typical outside of the @taylorize context. It also means the sequence of operations will be explicit for a compound expression rather than implicit.

  • Updating operators such as +=, *=, etc., are not supported. For example, the expression x += y is not recognized by @taylorize. Likewise, expressions such as x = x+y are not supported by @taylorize and should be replaced by equivalent expressions in which variables appear only on one side of the assignment; e.g. z = x+y; x = z. The introduction of such temporary variables z is left to the user.

  • The macro allows the use of array declarations through Array or Vector, but other ways (e.g. similar) are not yet implemented. Note that certain temporary arrays may be introduced to avoid re-computating certain expressions; only up-to-three indices expressions are currently handled.

  • Avoid using variables prefixed by an underscore, in particular _T, _S, _N and __idx, as well as ord; using them may lead to name collisions with some internal variables used in the constructed expressions.

  • Broadcasting is not recognized by @taylorize.

  • The macro may be used in combination with the common interface with DifferentialEquations.jl, for functions using the (du, u, p, t) in-place form, as we showed above. Other extensions allowed by DifferentialEquations may not be able to exploit it.

  • if-else blocks are recognized in their long form, but short-circuit conditional operators (&& and ||) are not. When comparing to a Taylor expansion, use operators such as iszero for if-else tests rather than comparing against numeric literals.

  • Input and output lengths should be determined at the time of @taylorize application (parse time), not at runtime. Avoid using the length of the input as an implicit indicator of whether to write all elements of the output. If conditional output of auxiliary equations is desired use explicit methods, such as through parameters or by setting auxiliary vector elements to zero, and assigning unwanted auxiliary outputs zero.

  • Expressions which correspond to function calls (so the head field is :call) which are not recognized by the parser are simply copied. The heuristics used, especially for vectors, may not work for all cases.

  • Use local for internal parameters, e.g., simple constant values; this improves performance. Do not use it if the variable is needed to be Taylor expanded during the integration step.

  • To examine the code generated for jetcoeffs! and _allocate_jetcoeffs! for a specific ODE function, follow the pendulum example above; create an expression by wrapping the ODE function (without @taylorize prefix) in a :()-block, and supply the expression to TaylorIntegration._make_parsed_jetcoeffs. This can help in debugging issues with either function generated by @taylorize.

  • @taylorize supports multi-threading via Threads.@threads. WARNING: this feature is experimental. Since thread-safety depends on the definition of each ODE, we cannot guarantee the resulting code to be thread-safe in advance. The user should check the resulting code to ensure that it is indeed thread-safe. For more information about multi-threading, the reader is referred to the Julia documentation.

We recommend to skim test/taylorize.jl, which implements different cases and highlights examples where the macro does not work, and how to solve the problem; read the information that is in the comments.

Please report any problems you may encounter.

+ end))

The first function has a similar structure as the hard-coded method of TaylorIntegration.jetcoeffs!, but uses low-level functions in TaylorSeries (e.g., sincos! above). Temporary Taylor1 objects as well as declared arrays are allocated once by TaylorIntegration._allocate_jetcoeffs!. More complex functions quickly become very difficult to read. Note that, if necessary, one can further optimize new_ex manually.

Limitations and some advice

The construction of the internal function obtained by using @taylorize is somewhat complicated and limited. Here we list some limitations and provide some advice.

  • Expressions must involve two arguments at most, which requires using parentheses appropriately: For example, res = a+b+c should be written as res = (a+b)+c. This may lead to more parentheses used in compound expressions than would be typical outside of the @taylorize context. It also means the sequence of operations will be explicit for a compound expression rather than implicit.

  • Updating operators such as +=, *=, etc., are not supported. For example, the expression x += y is not recognized by @taylorize. Likewise, expressions such as x = x+y are not supported by @taylorize and should be replaced by equivalent expressions in which variables appear only on one side of the assignment; e.g. z = x+y; x = z. The introduction of such temporary variables z is left to the user.

  • The macro allows the use of array declarations through Array or Vector, but other ways (e.g. similar) are not yet implemented. Note that certain temporary arrays may be introduced to avoid re-computating certain expressions; only up-to-three indices expressions are currently handled.

  • Avoid using variables prefixed by an underscore, in particular _T, _S, _N and __idx, as well as ord; using them may lead to name collisions with some internal variables used in the constructed expressions.

  • Broadcasting is not recognized by @taylorize.

  • The macro may be used in combination with the common interface with DifferentialEquations.jl, for functions using the (du, u, p, t) in-place form, as we showed above. Other extensions allowed by DifferentialEquations may not be able to exploit it.

  • if-else blocks are recognized in their long form, but short-circuit conditional operators (&& and ||) are not. When comparing to a Taylor expansion, use operators such as iszero for if-else tests rather than comparing against numeric literals.

  • Input and output lengths should be determined at the time of @taylorize application (parse time), not at runtime. Avoid using the length of the input as an implicit indicator of whether to write all elements of the output. If conditional output of auxiliary equations is desired use explicit methods, such as through parameters or by setting auxiliary vector elements to zero, and assigning unwanted auxiliary outputs zero.

  • Expressions which correspond to function calls (so the head field is :call) which are not recognized by the parser are simply copied. The heuristics used, especially for vectors, may not work for all cases.

  • Use local for internal parameters, e.g., simple constant values; this improves performance. Do not use it if the variable is needed to be Taylor expanded during the integration step.

  • To examine the code generated for jetcoeffs! and _allocate_jetcoeffs! for a specific ODE function, follow the pendulum example above; create an expression by wrapping the ODE function (without @taylorize prefix) in a :()-block, and supply the expression to TaylorIntegration._make_parsed_jetcoeffs. This can help in debugging issues with either function generated by @taylorize.

  • @taylorize supports multi-threading via Threads.@threads. WARNING: this feature is experimental. Since thread-safety depends on the definition of each ODE, we cannot guarantee the resulting code to be thread-safe in advance. The user should check the resulting code to ensure that it is indeed thread-safe. For more information about multi-threading, the reader is referred to the Julia documentation.

We recommend to skim test/taylorize.jl, which implements different cases and highlights examples where the macro does not work, and how to solve the problem; read the information that is in the comments.

Please report any problems you may encounter.