Library

Exported functions

TaylorIntegration.taylorinteg — Function.

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

General-purpose Taylor integrator for the explicit ODE $\dot{x}=f(x, p, t)$. The initial condition are specified by x0 at time t0, and any parameters encoded in params. The initial condition 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.

It returns a vector with the values of time (independent variable), and a vector with the computed values of the dependent variable(s). The integration stops when time is larger than tmax, in which case the last returned value(s) correspond to that time, 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 )

taylorinteg(f, x0, trange, order, abstol, params[=nothing]; keyword... )

General-purpose Taylor integrator for the explicit ODE $\dot{x}=f(t,x)$ with initial condition specified by x0::{T<:Number} or x0::Vector{T} at time t0.

The method returns a vector with the computed values of the dependent variable(s), evaluated only at the times specified by the range trange.

Examples

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

xv = taylorinteg(f!, [3, 3], 0.0:0.001:0.3, 25, 1.0e-20, maxsteps=100 );

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

taylorinteg(f, g, x0, trange, order, abstol; kwargs... )

Root-finding method of taylorinteg.

Given a function g(dx, x, params, t), called the event function, taylorinteg checks for the occurrence of a root of g evaluated at the solution; that is, it checks for the occurrence of an event or condition specified by g=0. Then, taylorinteg attempts to find that root (or event, or crossing) by performing a Newton-Raphson process. When called with the eventorder=n keyword argument, taylorinteg searches for the roots of the n-th derivative of g, which is computed via automatic differentiation.

The current keyword argument are:

maxsteps[=500]: maximum number of integration steps.
parse_eqs[=true]: use the specialized method of jetcoeffs! created with @taylorize.
eventorder[=0]: order of the derivative of g whose roots are computed.
newtoniter[=10]: maximum Newton-Raphson iterations per detected root.
nrabstol[=eps(T)]: allowed tolerance for the Newton-Raphson process; T is the common type of t0, tmax (or eltype(trange)) and abstol.

Examples:

using TaylorIntegration

function pendulum!(t, x, dx)
    dx[1] = x[2]
    dx[2] = -sin(x[1])
    nothing
end

g(dx, x, params, t) = x[2]

x0 = [1.3, 0.0]

# find the roots of `g` along the solution
tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20)

# find the roots of the 2nd derivative of `g` along the solution
tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, 0.0, 22.0, 28, 1.0E-20; eventorder=2)

# times at which the solution will be returned
tv = 0.0:1.0:22.0

# find the roots of `g` along the solution; return the solution *only* at each value of `tv`
xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, tv, 28, 1.0E-20)

# find the roots of the 2nd derivative of `g` along the solution; return the solution *only* at each value of `tv`
xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, tv, 28, 1.0E-20; eventorder=2)

TaylorIntegration.lyap_taylorinteg — Function.

lyap_taylorinteg(f!, q0, t0, tmax, order, abstol[, f!]; maxsteps::Int=500)

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.

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.

Internal

TaylorIntegration.__jetcoeffs! — Method.

__jetcoeffs!(::Val{bool}, f, t, x, params)
__jetcoeffs!(::Val{bool}, f, t, x, dx, xaux, params)

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

TaylorIntegration._capture_fn_args_body! — Function.

_capture_fn_args_body!(ex, vout::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.

TaylorIntegration._defs_preamble! — Function.

_defs_preamble!(preamble, fnargs, d_indx, v_newindx, v_arraydecl, v_preamb, d_decl, [inloop=false])

Returns a vector with expressions defining the auxiliary variables in the preamble; it may modify d_indx if new variables are introduced. v_preamb is for bookkeeping the introduced variables.

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

TaylorIntegration._make_parsed_jetcoeffs — Function.

_make_parsed_jetcoeffs( ex, debug=false )

This function constructs a new function, equivalent to the differential equations, which exploits the mutating functions of TaylorSeries.jl.

TaylorIntegration._newfnbody — Method.

_newfnbody(fnbody, fnargs, d_indx)

Returns a new (modified) body of the function, a priori unfolding the expression graph (AST) as unary and binary calls, and a vector (v_newindx) with the name of auxiliary indexed variables.

TaylorIntegration._newhead — Method.

_newhead(fn, fnargs)

Creates the head of the new method of jetcoeffs!. Here, fn is the name of the passed function and fnargs is a vector with its arguments (which are two or three).

TaylorIntegration._parse_newfnbody! — Function.

_parse_newfnbody!(ex, preex, v_vars, v_assign, d_indx, v_newindx, v_arraydecl, [inloop=false])

Parses ex (the new body of the function) replacing the expressions to use the mutating functions of TaylorSeries, and building the preamble preex. This is done by traversing recursively the args of ex, updating the bookkeeping vectors v_vars and v_assign. d_indx is used to substitute back the explicit indexed variables.

TaylorIntegration._preamble_body — Function.

_preamble_body(fnbody, fnargs, debug=false)

Returns the preamble, the body and a guessed return variable, which will be used to build the parsed function. fnbody is the expression with the body of the original function, fnargs is a vector of symbols of the original diferential equations function.

TaylorIntegration._recursionloop — Method.

_recursionloop(fnargs, retvar)

Build the expression for the recursion-loop.

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.

TaylorIntegration._replace_expr! — Method.

_replace_expr!(ex, preex, i, aalhs, aarhs, v_vars, d_indx, v_newindx)

Replaces the calls in ex.args[i], and updates preex with the definitions of the expressions, based on the the LHS (aalhs) and RHS (aarhs) of the base assignment. The bookkeeping vectors (v_vars, v_newindx) are updated. d_indx is used to bring back the indexed variables.

TaylorIntegration._replacecalls! — Method.

_replacecalls!(fnold, newvar, v_vars)

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

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.

TaylorIntegration._stepsize — Method.

_stepsize(aux1, epsilon, k)

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

TaylorIntegration.findroot! — Method.

findroot!(t, x, dx, g_val_old, g_val, 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_val_old is the last-before-current value of event function g; g_val is the current value of the event function g; 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.

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.

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.

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

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.

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

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.

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.

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; g_now represents the current value of event function g; eventorder is the order of the derivative of the event function g whose root we are trying to find. Returns true if g_old and g_now have different signs (i.e., if one is positive and the other one is negative). Otherwise, if g_old and g_now have the same sign or if either one of them are a nothing value, then returns false.

TaylorIntegration.taylorstep! — Method.

taylorstep!(f, t, x, t0, order, abstol, params, parse_eqs=true) -> δt
taylorstep!(f!, t, x, dx, xaux, t0, order, abstol, params, 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).

Index

TaylorIntegration.__jetcoeffs!
TaylorIntegration._capture_fn_args_body!
TaylorIntegration._defs_preamble!
TaylorIntegration._extract_parts
TaylorIntegration._make_parsed_jetcoeffs
TaylorIntegration._newfnbody
TaylorIntegration._newhead
TaylorIntegration._parse_newfnbody!
TaylorIntegration._preamble_body
TaylorIntegration._recursionloop
TaylorIntegration._rename_indexedvars
TaylorIntegration._replace_expr!
TaylorIntegration._replacecalls!
TaylorIntegration._second_stepsize
TaylorIntegration._stepsize
TaylorIntegration.findroot!
TaylorIntegration.jetcoeffs!
TaylorIntegration.jetcoeffs!
TaylorIntegration.lyap_jetcoeffs!
TaylorIntegration.lyap_taylorinteg
TaylorIntegration.lyap_taylorstep!
TaylorIntegration.nrconvergencecriterion
TaylorIntegration.stabilitymatrix!
TaylorIntegration.stepsize
TaylorIntegration.surfacecrossing
TaylorIntegration.taylorinteg
TaylorIntegration.taylorstep!