Library

taylorinteg(f, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... )
taylorinteg(f, x0, t0, tmax, order, abstol, Val(false), params[=nothing]; kwargs... )
taylorinteg(f, x0, t0, tmax, order, abstol, Val(true), params[=nothing]; kwargs... )

General-purpose Taylor integrator for the explicit ODE $\dot{x}=f(x, p, t)$, where p are the parameters encoded in params. The initial conditions are specified by x0 at time t0; x0 may be of type T<:Number or Vector{T}, with T including TaylorN{T}; the latter case is of interest for jet transport applications.

The equations of motion are specified by the function f; we follow the same convention of DifferentialEquations.jl to define this function, i.e., f(x, p, t) or f!(dx, x, p, t); see the examples below.

The functions return a vector with the values of time (independent variable), and a vector with the computed values of the dependent variable(s), and if the method used involves Val(true) it also outputs the Taylor polynomial solutions obtained at each time step. The integration stops when time is larger than tmax, in which case the last returned value(s) correspond to tmax, or when the number of saved steps is larger than maxsteps.

The integration method uses polynomial expansions on the independent variable of order order; the parameter abstol serves to define the time step using the last two Taylor coefficients of the expansions. Make sure you use a large enough order to assure convergence.

Currently, the recognized keyword arguments are:

• maxsteps[=500]: maximum number of integration steps.
• parse_eqs[=true]: use the specialized method of jetcoeffs! created with @taylorize.

Examples

For one dependent variable the function f defines the RHS of the equation of motion, returning the value of $\dot{x}$. The arguments of this function are (x, p, t), where x are the dependent variables, p are the paremeters and t is the independent variable.

For several (two or more) dependent variables, the function f! defines the RHS of the equations of motion, mutating (in-place) the (preallocated) vector with components of $\dot{x}$. The arguments of this function are (dx, x, p, t), where dx is the preallocated vector of $\dot{x}$, x are the dependent variables, p are the paremeters entering the ODEs and t is the independent variable. The function may return this vector or simply nothing.

using TaylorIntegration

f(x, p, t) = x^2

tv, xv = taylorinteg(f, 3, 0.0, 0.3, 25, 1.0e-20, maxsteps=100 )

function f!(dx, x, p, t)
    for i in eachindex(x)
        dx[i] = x[i]^2
    end
    return nothing
end

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) )

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.

@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.

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

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}}}

__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).

_allocated_defs!(new_jetcoeffs, bkkeep)

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

_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.

_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.

_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.

_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).

_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.

_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.

_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).

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.

_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.

_recursionloop(fnargs, bkkeep)

Build the expression for the recursion-loop.

_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.

_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.

_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).

_returned_expr(bkkeep)

Constructs the expression to be returned by TaylorIntegration._allocate_jetcoeffs!

_second_stepsize(x, epsilon)

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

_split_arraydecl!(bkkeep)

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

_stepsize(aux1, epsilon, k)

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

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.

inbookkeeping(v, bkkeep::BookKeeping)

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

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.

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.

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!.

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.

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.

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.

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.

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.

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).