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:
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\),
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:
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:
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):
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 usesymengine.sin
– in contrast tomath.sin
ornumpy.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.) Seelarge_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='140091931032936'>¶ 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='140091931045728'>¶ 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 thei
-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 bey(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 iff_sym
is a generator function andn
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 givef_sym
as an argument, but in this case you must given
. 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. IfNone
, this will be automatically disabled forn>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 off
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
orsolve_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
andgenerate_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 orFalse
, 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
orLambdifyCSE
.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. IfNone
, this will be automatically disabled forn>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
orLambdifyCSE
. If the symbolic Jacobian has not been generated, it is generated by callinggenerate_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.
-
y_dict
¶ The current state of the system as a dictionary mapping dynamical variables to their current value. Note that if you use this often, you may want to use
self.y
instead for efficiency.
-
set_initial_value
(initial_value, time=0.0)¶ Same as the analogous function in SciPy’s ODE, except that it also accepts the initial_value in form of a dictionary that maps dynamical variables to their initial value.
-
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 viaode
"RK45"
– Dormand’s and Prince’s explicit fifth-order method viasolve_ivp
"dop853"
– DoP853 (explicit) viaode
"RK23"
– Bogacki’s and Shampine’s explicit third-order method viasolve_ivp
"BDF"
– Implicit backward-differentiation formula viasolve_ivp
"lsoda"
– LSODA (implicit) viaode
"LSODA"
– LSODA (implicit) viasolve_ivp
"Radau"
– The implicit Radau method viasolve_ivp
"vode"
– VODE (implicit) viaode
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
andset_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 dimensionn
- 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.
- negative arguments of
-
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 usingcompile_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). IfNone
, this will be automatically disabled forn>10
.
-
set_initial_value
(y, time=0.0)¶ Same as the analogous function in SciPy’s ODE, except that it also accepts the initial_value in form of a dictionary that maps dynamical variables to their initial value.
-
set_integrator
(name, interpolate=None, **kwargs)¶ Same as for
jitcode
except thatinterpolate
defaults toFalse
for RK45 and Radau to avoid an accumulation error through the nature of the interpolation algorithm. UsingLSODA
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
’sintegrate
andode
’sintegrate
.- 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.
-
y_dict
¶ The current state of the system as a dictionary mapping dynamical variables to their current value. Note that if you use this often, you may want to use
self.y
instead for efficiency.
-
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). IfNone
, this will be automatically disabled forn>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 thegroups
argument of the constructor).
-
set_integrator
(name, interpolate=False, **kwargs)¶ Same as for
jitcode
except thatinterpolate
defaults toFalse
for RK45 and Radau to avoid an accumulation error through the nature of the interpolation algorithm. UsingLSODA
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]
-
y_dict
¶ The current state of the system as a dictionary mapping dynamical variables to their current value. Note that if you use this often, you may want to use
self.y
instead for efficiency.
-
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)
|