The JiTCODE module

Note and remember that some relevant information can be found in the common JiTC*DE documentation. This includes installation instructions, compiler issues and optimisation, general performance considerations, how to implement network dynamics, and a small FAQ.

Overview

JiTCODE (just-in-time compilation for ordinary differential equations) is an extension of SciPy’s ODE (scipy.integrate.ode) or Solve IVP (scipy.integrate.solve_ivp). Where the latter take a Python function as an argument, JiTCODE takes an iterable (or generator function or dictionary) of symbolic expressions, which it translates to C code, compiles on the fly, and uses as the function to feed into SciPy’s ODE or Solve IVP. Symbolic expressions are mostly handled by SymEngine, SymPy’s compiled-backend-to-be (see SymPy vs. SymEngine for details).

This design has the following advantages:

  • Speed boost through compilation There are two main computational burdens when numerically solving ODEs. The first are simple vector operations for the integration procedure, which can be efficiently handled by NumPy or similar. The second is calculating the derivative (and Jacobian, if desired), which can be efficiently handled with just-in-time compilation. Note that for smaller ODEs and small integration times, the Python overhead may still take a considerable toll in comparison to pure compiled code.
  • Speed boost through symbolic optimisation If your derivative is automatically generated by some routine, you can simplify it symbolically to boost the speed. In fact, blatant optimisations such as \(y·(x-x)=0\) are done on the fly by SymEngine. This is for example interesting if you want to simulate dynamics on a sparse network, as non-existing links are not taken into account when calculating the derivative when integrating.
  • Automatically calculated Jacobian and Lyapunov exponents As the derivative is provided symbolically, SymEngines’s automatic differentiation routines can be employed to calculate the Jacobian desired by some integrators as well as the ODE for the tangent vector required for calculating the Lyapunov exponents (see Calculating Lyapunov exponents with jitcode_lyap).
  • Symbolic interface You can enter your differential equations almost like you would on paper. Also, if you are working with SymPy or SymEngine anyway – e.g., to calculate fixed points –, you do not need to bother much with translating your equations.

If compilation fails to work for whatever reason, pure Python functions can be employed as a fallback.

For most applications, the only difference to SciPy’s ODE in terms of handling is that you have to supply the derivative in the correct format – JiTCODE automatically takes care of the compilation. However, you can also tweak this step by step, if you desire (see Details of the building process and the notes on compilation).

As with SciPy’s ODE, this documentation assumes that the differential equation you want to solve is:

\[\dot{y} = f(t,y)\]

A quick example

This example is particularly intended for readers familiar with SciPy’s ODE. If you aren’t, skip to the next section.

Suppose our differential equation is \(\dot{y} = f(y)\) with \(y∈ℝ^4\),

\[\begin{split}f(y) = \left( \begin{matrix} y_0 ( a-y_0 ) ( y_0-1) - y_1 + k (y_2 - y_0) \\ b_1 y_0 - c y_1 \\ y_2 ( a-y_2 ) ( y_2-1 ) - y_3 + k (y_0 - y_2)\\ b_2 y_2 - c y_3 \end{matrix} \right),\end{split}\]

and \(a = -0.025794\), \(b_1 = 0.0065\), \(b_2 = 0.0135\), \(c = 0.02\), and \(k = 0.128\). Then the following code integrates the above for 100000 time units (after discarding 2000 time units of transients), with \(y(t=0) = (1,2,3,4)\), and writes the results to timeseries.dat:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
from jitcode import jitcode, y
import numpy as np

a  = -0.025794
b1 =  0.0065
b2 =  0.0135
c  =  0.02
k  =  0.128

f = [
	y(0) * ( a-y(0) ) * ( y(0)-1.0 ) - y(1) + k * (y(2) - y(0)),
	b1*y(0) - c*y(1),
	y(2) * ( a-y(2) ) * ( y(2)-1.0 ) - y(3) + k * (y(0) - y(2)),
	b2*y(2) - c*y(3)
	]

initial_state = np.array([1.,2.,3.,4.])

ODE = jitcode(f)
ODE.set_integrator("dopri5")
ODE.set_initial_value(initial_state,0.0)

times = 2000+np.arange(100000)
data = []
for time in times:
	data.append(ODE.integrate(time))

np.savetxt("timeseries.dat", data)

An explicit example

This example is particularly intended for readers who are using a numerical integrator for the first time.

Suppose, we want to implement the Lotka–Volterra model, which is described by the following equations:

\[\begin{split}\begin{alignat*}{3} \dot{B} &=& γ · B &- φ · R · B\\ \dot{R} &=&\, -ω · R &+ ν · R · B \end{alignat*}\end{split}\]

with \(γ = 0.6\), \(φ = 1.0\), \(ω = 0.5\), and \(ν = 0.5\).

We start with a few imports:

from jitcode import y, jitcode
import numpy as np

… and defining the cotrol parameters:

γ = 0.6
φ = 1.0
ω = 0.5
ν = 0.5

The y that we imported from jitcode has to be used for the dynamical variables. However, to make our code use the same notation as the above equation, we can rename them:

R,B = [ y(i) for i in range(2) ]

We here might as well have written R,B = y(0),y(1), but the above is more handy for larger systems. Note that the following code is written such that the order of our dynamical variables only matters for the output.

We implement the differential equation as a dictionary, mapping each of the dynamical variables to its derivative:

lotka_volterra_diff = {
		B:  γ*B - φ*R*B,
		R: -ω*R + ν*R*B,
	}

Note that there are other ways to define the derivative, e.g., used in the previous and following example.

Now, we have everything, we need to instantiate jitcode, i.e., to initialise the integration object:

ODE = jitcode(lotka_volterra_diff)

Before we start integrating, we have to choose an integrator and define initial conditions. We here choose the 5th-order Dormand–Prince method, random initial states (between 0 and 1) and start at \(t=0\).

ODE.set_integrator("dopri5")
initial_state = np.array(np.random.random(2))
ODE.set_initial_value(initial_state,0.0)

We then define an array of time points where we want to observe the system and an empty list that shall be filled with our results.

times = np.arange(0.0,100,0.1)
data = []

Finally, we perform the actual integration by calling ode.integrate for each of our times. This returns the state of the system after integration, which we put into a tuple together with the current time and append to our data list. (The asterisk (*) unpacks an iterable. We might as well have written [time,state[0],state[1]] instead of [time,*state].)

for time in times:
	state = ODE.integrate(time)
	data.append( [time,*state] )

We can now plot or analyse our data, but that’s beyond the scope of JiTCODE. So we just save it:

np.savetxt("timeseries.dat", data)

Taking everything together, our code is:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
from jitcode import y, jitcode
import numpy as np

γ = 0.6
φ = 1.0
ω = 0.5
ν = 0.5

R,B = [ y(i) for i in range(2) ]

lotka_volterra_diff = {
		B:  γ*B - φ*R*B,
		R: -ω*R + ν*R*B,
	}

ODE = jitcode(lotka_volterra_diff)
ODE.set_integrator("dopri5")
initial_state = np.array(np.random.random(2))
ODE.set_initial_value(initial_state,0.0)

times = np.arange(0.0,100,0.1)
data = []
for time in times:
	state = ODE.integrate(time)
	data.append( [time,*state] )

np.savetxt("timeseries.dat", data)

Details of the building process

To generate the functions needed by SciPy’s ODE or Solve IVP, JiTCODE executes a series of distinct processing steps, each of which can be invoked through a command and tweaked as needed. These commands execute previous steps as needed, i.e., if they have not been performed yet. If you are happy with the default options, however, you do not need to bother with this and can just use the commands at the very end of the chain, namely set_integrator, set_initial_value, or save_compiled.

The following diagram details which command calls which other command when needed:

digraph foo { bgcolor = "transparent" edge [arrowsize=0.5]; node [shape=box fontname="Monospace" fontsize=10 margin=0.05 fillcolor="#C8D5E3" penwidth=0 style=filled]; {rank = same; "save_compiled"; "set_integrator"} "set_initial_value" -> "set_integrator"; "save_compiled" -> "compile_C"; "set_integrator" -> "generate_functions" -> "compile_C" -> "generate_f_C" -> "generate_helpers_C"; "compile_C" -> "generate_jac_C" -> "generate_jac_sym"; "generate_jac_C" -> "generate_helpers_C"; "generate_functions" -> "generate_lambdas" -> "generate_f_lambda"; "generate_lambdas" -> "generate_jac_lambda"-> "generate_jac_sym"; }

A more complicated example

This example showcases several advanced features of JiTCODE that are relevant for an efficient integration of more complex systems as well as how to deal with some special situations. Therefore it is pretty bizarre from a dynamical perspective.

Suppose we want to integrate a system of \(N=500\) Rössler oscillators, with the \(i\)-th oscillator being described by the following differential equations (note that we used \(v\) instead of \(y\) for the second dynamical variable to avoid overloading symbols):

\[\begin{split}\begin{alignedat}{1} \dot{x}_i &= -ω_i v_i - z_i + k \sum_{j=0}^N A_{ij} (x_j-x_i) \\ \dot{v}_i &= ω_i x_i + a v_i \\ \dot{z}_i &= b + z_i (x_i -c) + k \sin(t) \sum_{j=0}^N (x_j-x_i) \end{alignedat}\end{split}\]

The control parameters shall be \(a = 0.165\), \(b = 0.2\), \(c = 10.0\), and \(k = 0.01\). The (frequency) parameter \(ω_i\) shall be picked randomly from the uniform distribution on \([0.8,1.0]\) for each \(i\). \(A∈ℝ^{N×N}\) shall be the adjacency matrix of a one-dimensional small-world network (which shall be provided by a function small_world_network in the following example code). So, the \(x\) compenents are coupled diffusively with a small-world coupling topology, while the \(z\) components are coupled diffusively to their mean field, with the coupling term being modulated with \(\sin(t)\).

Without further ado, here is the example code (complete running example); highlighted lines will be commented upon below:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
from jitcode import jitcode, y, t
import numpy as np
import symengine

# parameters
# ----------

N = 500
ω = np.random.uniform(0.8,1.0,N)
a = 0.165
b = 0.2
c = 10.0
k = 0.01

# get adjacency matrix of a small-world network
A = small_world_network(
	number_of_nodes = N,
	nearest_neighbours = 20,
	rewiring_probability = 0.1
	)

# generate differential equations
# -------------------------------

sum_z = symengine.Symbol("sum_z")
helpers = [( sum_z, sum( y(3*j+2) for j in range(N) ) )]

def f():
	for i in range(N):
		coupling_sum = sum( y(3*j)-y(3*i) for j in range(N) if A[i,j] )
		coupling_term = k * symengine.sin(t) * coupling_sum
		yield -ω[i] * y(3*i+1) - y(3*i+2) + coupling_term
		yield  ω[i] * y(3*i) + a*y(3*i+1)
		coupling_term_2 = k * (sum_z-N*y(3*i+2))
		yield b + y(3*i+2) * (y(3*i) - c) + coupling_term_2

# integrate
# ---------

initial_state = np.random.random(3*N)

ODE = jitcode(f, helpers=helpers, n=3*N)
ODE.generate_f_C(simplify=False, do_cse=False, chunk_size=150)
ODE.set_integrator('dopri5')
ODE.set_initial_value(initial_state,0.0)

# data structure: x[0], v[0], z[0], x[1], …, x[N], v[N], z[N]
data = np.vstack(ODE.integrate(T) for T in range(10,100000,10))

Explanation of selected features and choices:

  • The values of \(ω\) are initialised globally (line 9). We cannot just define a function here, because the parameter is used twice for each oscillator. Moreover, if we were trying to calculate Lyapunov exponents or the Jacobian, the generator function would be called multiple times, and thus the value of the parameter would not be consistent (which would be desastrous).
  • Since we need \(\sum_{j=0}^N x_j\) to calculate the derivative of \(z\) for every oscillator, it is prudent to only calculate this once. Therefore we define a helper symbol for this in lines 25 and 26, which we employ in line 34. (See the arguments of jitcode for details.)
  • In line 31, we implement \(\sin(t)\). For this purpose we had to import t in line 1. Also, we need to use symengine.sin – in contrast to math.sin or numpy.sin.
  • As this is a large system, we use a generator function instead of a list to define \(f\) (lines 28-35) and have the code automatically be split into chunks of 150 lines, corresponding to the equations of fifty oscillators, in line 43. (For this system, any reasonably sized multiple of 3 is a good choice for chunk_size; for other systems, the precise choice of the value may be more relevant.) See large_systems for the rationale.

Calculating Lyapunov exponents with jitcode_lyap

jitcode_lyap is a simple extension of jitcode that almost automatically handles calculating Lyapunov exponents by evolving tangent vectors [BGGS80]. It works just like jitcode, except that it generates and integrates additional differential equations for the tangent vectors. After every call of integrate, the tangent vectors are orthonormalised, and the “local” Lyapunov exponents for this integration step are returned alongside with the system’s state. These can then be further processed to obtain the Lyapunov exponents. The tangent vectors are initialised with random vectors, and you have to take care of the preiterations that the tangent vectors require to align themselves.

Estimates for the Lyapunov vectors are returned as well.

For instance, we can calculate and print the Lyapunov exponents for the system from A quick example as follows (changes highlighted):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
from jitcode import jitcode_lyap, y
from scipy.stats import sem
import numpy as np

a  = -0.025794
b1 =  0.0065
b2 =  0.0135
c  =  0.02
k  =  0.128

f = [
	y(0) * ( a-y(0) ) * ( y(0)-1.0 ) - y(1) + k * (y(2) - y(0)),
	b1*y(0) - c*y(1),
	y(2) * ( a-y(2) ) * ( y(2)-1.0 ) - y(3) + k * (y(0) - y(2)),
	b2*y(2) - c*y(3)
	]

initial_state = np.array([1.,2.,3.,4.])

n = len(f)
ODE = jitcode_lyap(f, n_lyap=n)
ODE.set_integrator("vode")
ODE.set_initial_value(initial_state,0.0)

times = range(10,100000,10)
lyaps = []
for time in times:
	lyaps.append(ODE.integrate(time)[1])

# converting to NumPy array for easier handling
lyaps = np.vstack(lyaps)

for i in range(n):
	lyap = np.average(lyaps[1000:,i])
	stderr = sem(lyaps[1000:,i]) # Note that this only an estimate
	print("%i. Lyapunov exponent: % .4f ± %.4f" % (i+1,lyap,stderr))

Calculating transversal Lyapunov exponents with jitcode_transversal_lyap

jitcode_transversal_lyap is a variant of jitcode_lyap that calculates Lyapunov exponents transversal to a user-defined synchronisation manifold. It automatically conflates the differential equations of a group of synchronised components into one equation on the synchronisation manifold. Moreover, it transforms the equations for the tangent vectors such that the tangent vector is automatically orthogonal to the synchronisation manifold. If you are interested in the mathematical details, please read the accompanying paper.

For instance, let’s interpret the system from A quick example as two oscillators (which is what it is), one consisting of the first and second and one of the third and fourth component. Furthermore, let’s change the control parameters a bit to make the two oscillators identical. We can then calculate the transversal Lyapunov exponents to the synchronisation manifold as follows (important changes are highlighted):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
from jitcode import jitcode_transversal_lyap, y
from scipy.stats import sem
import numpy as np

a = -0.025794
b =  0.01
c =  0.02
k =  0.128

f = [
	y(0) * ( a-y(0) ) * ( y(0)-1.0 ) - y(1) + k * (y(2) - y(0)),
	b*y(0) - c*y(1),
	y(2) * ( a-y(2) ) * ( y(2)-1.0 ) - y(3) + k * (y(0) - y(2)),
	b*y(2) - c*y(3)
	]

initial_state = np.array([1.,2.])

groups = [ [0,2], [1,3] ]
ODE = jitcode_transversal_lyap(f, groups=groups)
ODE.set_integrator("lsoda")
ODE.set_initial_value(initial_state,0.0)

times = range(10,100000,10)
lyaps = []
for time in times:
	lyaps.append(ODE.integrate(time)[1])

lyap = np.average(lyaps[1000:])
stderr = sem(lyaps[1000:]) # Note that this only an estimate
print("transversal Lyapunov exponent: % .4f ± %.4f" % (lyap,stderr))

Note that the initial state (line 17) is reduced in dimensionality and there is only one component for each synchronisation group.

Command reference

y = <MagicMock name='mock.Function()' id='139780554998672'>

the symbol for the state that must be used to define the differential equation. It is a function and the integer argument denotes the component. You may just as well define an analogous function directly with SymEngine or SymPy, but using this function is the best way to get the most of future versions of JiTCODE, in particular avoiding incompatibilities.

t = <MagicMock name='mock.Symbol()' id='139780555015560'>

the symbol for time for defining the differential equation. If your differential equation has no explicit time dependency (“autonomous system”), you do not need this. You may just as well define an analogous symbol directly with SymEngine or SymPy, but using this function is the best way to get the most of future versions of JiTCODE, in particular avoiding incompatibilities.

class UnsuccessfulIntegration

This exception is raised when the integrator cannot meet the accuracy and step-size requirements. If you want to know the exact state of your system before the integration fails or similar, catch this exception.

The main class

class jitcode(f_sym=(), helpers=None, wants_jacobian=False, n=None, control_pars=(), verbose=True, module_location=None)
Parameters:
f_sym : iterable of symbolic expressions or generator function yielding symbolic expressions or dictionary

If an iterable or generator function, the i-th element is the i-th component of the value of the ODE’s derivative \(f(t,y)\). If a dictionary, it has to map the dynamical variables to its derivatives and the dynamical variables must be y(0), y(1), .

helpers : list of length-two iterables, each containing a symbol and a symbolic expression

Each helper is a variable that will be calculated before evaluating the derivative and can be used in the latter’s computation. The first component of the tuple is the helper’s symbol as referenced in the derivative or other helpers, the second component describes how to compute it from t, y and other helpers. This is for example useful to realise a mean-field coupling, where the helper could look like (mean, sum(y(i) for i in range(100))/100). (See A more complicated example for an example.)

wants_jacobian : boolean

Tell JiTCODE to calculate and compile the Jacobian. For vanilla use, you do not need to bother about this as this is automatically set to True if the selected method of integration desires the Jacobian. However, it is sometimes useful if you want to manually apply some code-generation steps (e.g., to apply some tweaks).

n : integer

Length of f_sym. While JiTCODE can easily determine this itself (and will, if necessary), this may take some time if f_sym is a generator function and n is large. Take care that this value is correct – if it isn’t, you will not get a helpful error message.

control_pars : iterable of symbols

Each symbol corresponds to a control parameter that can be used when defining the equations and set after compilation using set_parameters (in the same order as given here). Using this makes sense if you need to do a parameter scan with short integrations for each parameter and you spend a considerable amount of time compiling.

verbose : boolean

Whether JiTCODE shall give progress reports on the processing steps.

module_location : string

location of a module file from which functions are to be loaded (see save_compiled). If you use this, you need not give f_sym as an argument, but in this case you must give n. Depending on the arguments you provide, functionalities such as recompiling may not be available; but then the entire point of this option is to avoid these.

generate_jac_sym(simplify=True)

generates the Jacobian using SymEngine’s differentiation.

Parameters:
simplify : boolean

Whether the resulting Jacobian should be simplified (with ratio=1.0). This is almost always a good thing.

generate_f_C(simplify=None, do_cse=False, chunk_size=100)

translates the derivative to C code using SymEngine’s C-code printer.

Parameters:
simplify : boolean or None

Whether the derivative should be simplified (with ratio=1.0) before translating to C code. The main reason why you could want to enable this is if you expect your derivative not to be optimised and not be so large that simplifying takes a considerable amount of time. If None, this will be automatically disabled for n>10.

do_cse : boolean

Whether SymPy’s common-subexpression detection should be applied before translating to C code. It is almost always better to let the compiler do this (unless you want to set the compiler optimisation to -O2 or lower): For simple differential equations this should not make any difference to the compiler’s optimisations. For large ones, it may make a difference but also take long. As this requires all entries of f at once, it may void advantages gained from using generator functions as an input. Also, this feature uses SymPy and not SymEngine.

chunk_size : integer

If the number of instructions in the final C code exceeds this number, it will be split into chunks of this size. See Handling very large differential equations on why this is useful and how to best choose this value. If smaller than 1, no chunking will happen.

generate_jac_C(do_cse=False, chunk_size=100, sparse=True)

translates the symbolic Jacobian to C code using SymEngine’s C-code printer. If the symbolic Jacobian has not been generated, it generates it by calling generate_jac_sym.

Parameters:
do_cse : boolean

Whether SymPy’s common-subexpression detection should be applied before translating to C code. It is almost always better to let the compiler do this (unless you want to set the compiler optimisation to -O2 or lower): For simple differential equations this should not make any difference to the compiler’s optimisations. For large ones, it may make a difference but also take long. As this requires the entire Jacobian at once, it may void advantages gained from using generator functions as an input. Also, this feature uses SymPy and not SymEngine.

chunk_size : integer

If the number of instructions in the final C code exceeds this number, it will be split into chunks of this size. See Handling very large differential equations on why this is useful and how to best choose this value. If smaller than 1, no chunking will happen.

sparse : boolean

Whether a sparse Jacobian should be assumed for optimisation. Note that this does not mean that the Jacobian is stored, parsed or handled as a sparse matrix. This kind of optimisation would require ode or solve_ivp to be able to handle sparse matrices without structure in the sparseness.

generate_helpers_C(chunk_size=100)

translates the helpers to C code using SymEngine’s C-code printer.

Parameters:
chunk_size : integer

If the number of instructions in the final C code exceeds this number, it will be split into chunks of this size. See large_systems on why this is useful.

If there is an obvious grouping of your helpers, the group size suggests itself for chunk_size.

If smaller than 1, no chunking will happen.

compile_C(extra_compile_args=None, extra_link_args=None, verbose=False, modulename=None, omp=False)

compiles the C code (using Setuptools) and loads the compiled functions. If no C code exists, it is generated by calling generate_f_C and generate_jac_C. For detailed information on the arguments and other ways to tweak the compilation, read these notes.

Parameters:
extra_compile_args : iterable of strings
extra_link_args : iterable of strings

Arguments to be handed to the C compiler or linker, respectively.

verbose : boolean

Whether the compiler commands shall be shown. This is the same as Setuptools’ verbose setting.

modulename : string or None

The name used for the compiled module.

omp : pair of iterables of strings or boolean

What compiler arguments shall be used for multiprocessing (using OpenMP). If True, they will be selected automatically. If empty or False, no compilation for multiprocessing will happen (unless you supply the relevant compiler arguments otherwise).

generate_f_lambda(simplify=None, do_cse=False)

translates the symbolic derivative to a function using SymEngines Lambdify or LambdifyCSE.

Parameters:
simplify : boolean

Whether the derivative should be simplified (with ratio=1.0) before translating to C code. The main reason why you could want to enable this is if you expect your derivative not to be optimised and not be so large that simplifying takes a considerable amount of time. If None, this will be automatically disabled for n>10.

do_cse : boolean

Whether a common-subexpression detection, namely LambdifyCSE, should be used.

generate_jac_lambda(do_cse=False)

translates the symbolic Jacobain to a function using SymEngines Lambdify or LambdifyCSE. If the symbolic Jacobian has not been generated, it is generated by calling generate_jac_sym.

Parameters:
do_cse : boolean

Whether a common-subexpression detection, namely LambdifyCSE, should be used.

generate_lambdas()

If they do not already exists, this generates lambdified functions by calling self.generate_f_lambda() and, if wanted, generate_jac_lambda().

generate_functions()

The central function-generating function. Tries to compile the derivative and, if wanted, the Jacobian. If this fails, it generates lambdified functions as a fallback.

set_initial_value(initial_value, time=0.0)

Same as the analogous function in SciPy’s ODE.

set_integrator(name, nsteps=1000000, interpolate=True, **integrator_params)

Analogous to the function in SciPy’s ODE with the same name. This automatically generates the derivative and Jacobian if they do not exist yet and are needed. You can also choose integrators from scipy.integrate.solve_ivp.

Parameters:
name: name of the integrator

One of the following (or a new method supported by either backend):

  • "dopri5" – Dormand’s and Prince’s explicit fifth-order method via ode
  • "RK45" – Dormand’s and Prince’s explicit fifth-order method via solve_ivp
  • "dop853" – DoP853 (explicit) via ode
  • "RK23" – Bogacki’s and Shampine’s explicit third-order method via solve_ivp
  • "BDF" – Implicit backward-differentiation formula via solve_ivp
  • "lsoda" – LSODA (implicit) via ode
  • "LSODA" – LSODA (implicit) via solve_ivp
  • "Radau" – The implicit Radau method via solve_ivp
  • "vode" – VODE (implicit) via ode

The solve_ivp methods are usually slightly faster for large differential equations, but they come with a massive overhead that makes them considerably slower for small differential equations. Implicit solvers are slower than explicit ones, except for stiff problems. If you don’t know what to choose, start with "dopri5".

nsteps: integer

Same as the respective parameter of the ode solvers, but with a higher default value to avoid annoying errors when getting rid of transients.

interpolate: boolean

Whether the sampled solutions for solve_ivp solvers shall be obtained using interpolation. If your sampling step is small, this may make things faster; otherwise it depends. This may also make the results slightly less accurate.

integrator_params

Parameters passed to the respective integrator. See its documentation for more.

set_parameters(*args)

Same as set_f_params and set_jac_params for SciPy’s ODE (both sets of parameters are set simultaneuosly, because they should be the same anyway).

check(fail_fast=True)

Performs a series of checks that may not be feasible at runtime (usually due to their length). Whenever you run into an error that you cannot make sense of, try running this. It checks for the following mistakes:

  • negative arguments of y
  • arguments of y that are higher than the system’s dimension n
  • unused variables
Parameters:
fail_fast : boolean

whether to abort on the first failure. If false, an error is raised only after all problems are printed.

render_and_write_code(expressions, name, chunk_size=100, arguments=(), omp=True)

Writes expressions to code.

Parameters:
expressions: iterator

expressions to be written

name: string

unique name of what is computed

chunk_size: integer

size of chunks. If smaller than 1, no chunking happens.

arguments: list of tuples

Each tuple contains the name, type, and size (optional, for arrays) of an argument needed by the code. This is so the arguments can be passed to chunked functions.

omp: boolean

whether OMP pragmas should be included

save_compiled(destination='', overwrite=False)

saves the module file with the compiled functions for later use (see the module_location argument). If no compiled derivative exists, it tries to compile it first using compile_C. In most circumstances, you should not rename this file, as the filename is needed to determine the module name.

Parameters:
destination : string specifying a path

If this specifies only a directory (don’t forget the trailing / or similar), the module will be saved to that directory. If empty (default), the module will be saved to the current working directory. Otherwise, the functions will be (re)compiled to match that filename. A file ending will be appended if needed.

overwrite : boolean

Whether to overwrite the specified target if it already exists.

Returns:
filename : string

The destination that was actually used.

Lyapunov exponents

class jitcode_lyap(f_sym=(), n_lyap=-1, simplify=None, **kwargs)

the handling is the same as that for jitcode except for:

Parameters:
n_lyap : integer

Number of Lyapunov exponents to calculate. If negative or larger than the dimension of the system, all Lyapunov exponents are calculated.

simplify : boolean

Whether the differential equations for the tangent vector shall be subjected to SymEngine’s simplify. Doing so may speed up the time evolution but may slow down the generation of the code (considerably for large differential equations). If None, this will be automatically disabled for n>10.

set_initial_value(y, time=0.0)

Same as the analogous function in SciPy’s ODE.

set_integrator(name, interpolate=None, **kwargs)

Same as for jitcode except that interpolate defaults to False for RK45 and Radau to avoid an accumulation error through the nature of the interpolation algorithm. Using LSODA as an integrator is discouraged.

integrate(*args, **kwargs)

Like SciPy’s ODE’s integrate, except for orthonormalising the tangent vectors and:

Returns:
y : one-dimensional NumPy array

The state of the system. Same as the output of jitcode’s integrate and ode’s integrate.

lyaps : one-dimensional NumPy array

The “local” Lyapunov exponents as estimated from the growth or shrinking of the tangent vectors during the integration time of this very integrate command, i.e., \(\frac{\ln (α_i^{(p)})}{s_i}\) in the notation of [BGGS80]

lyap_vectors : list of one-dimensional NumPy arrays

The Lyapunov vectors (normalised) after integration.

class jitcode_transversal_lyap(f_sym=(), groups=(), simplify=None, **kwargs)

Calculates the largest Lyapunov exponent in orthogonal direction to a predefined synchronisation manifold, i.e. the projection of the tangent vector onto that manifold vanishes. In contrast to jitcode_restricted_lyap, this performs some transformations tailored to this specific application that may strongly reduce the number of differential equations and ensure a dynamics on the synchronisation manifold.

The handling is the same as that for jitcode except for:

Parameters:
groups : iterable of iterables of integers

each group is an iterable of indices that identify dynamical variables that are synchronised on the synchronisation manifold.

simplify : boolean

Whether the transformed differential equations shall be subjected to SymEngine’s simplify. Doing so may speed up the time evolution but may slow down the generation of the code (considerably for large differential equations). If None, this will be automatically disabled for n>10.

set_initial_value(y, time=0.0)

Like the analogous function of jitcode/scipy.integrate.ode, except that only one initial value per group of synchronised components has to be provided (in the same order as the groups argument of the constructor).

set_integrator(name, interpolate=False, **kwargs)

Same as for jitcode except that interpolate defaults to False for RK45 and Radau to avoid an accumulation error through the nature of the interpolation algorithm. Using LSODA as an integrator is discouraged.

integrate(*args, **kwargs)

Like SciPy’s ODE’s integrate, except for orthonormalising the tangent vector and:

Returns:
y : one-dimensional NumPy array

The state of the system. Only one initial value per group of synchronised components is returned (in the same order as the groups argument of the constructor).

lyaps : one-dimensional NumPy array

The “local” Lyapunov exponent as estimated from the growth or shrinking of the tangent vector during the integration time of this very integrate command, i.e., \(\frac{\ln (α_i^{(p)})}{s_i}\) in the notation of [BGGS80]

class jitcode_restricted_lyap(f_sym=(), vectors=(), **kwargs)

Calculates the largest Lyapunov exponent in orthogonal direction to a predefined plane, i.e. the projection of the tangent vector onto that plane vanishes. See this test for an example of usage. The handling is the same as that for jitcode_lyap except for:

Parameters:
vectors : iterable of pairs of NumPy arrays

A basis of the plane, whose projection shall be removed.

What doesn’t work (yet)

The following feature of SciPy’s ODE or Solve IVP cannot be used through JiTCODE:

  • The zvode integrator and complex arithmetics in general, as they are not easy implementable in C.

References

[BGGS80](1, 2, 3)
  1. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn: Lyapunov Characteristic Exponents for smooth dynamical systems and for Hamiltonian systems; A method for computing all of them. Meccanica 15, pp. 9–30 (1980), 10.1007/BF02128236.