Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@

# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = sphinx-build
SPHINXBUILD = python $(shell which sphinx-build)
PAPER =
BUILDDIR = _build

Expand Down
15 changes: 7 additions & 8 deletions examples/kuramoto_network.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,21 +15,21 @@
:linenos:
:start-after: example-st\u0061rt
:dedent: 1
:emphasize-lines: 9, 10, 15, 20, 21, 24
:emphasize-lines: 11, 12, 17, 22, 23, 26

Explanation of selected features and choices:

* Line 9 is just a quick way to generate the network described above. For more complex networks, you will either have to write more complex function or use dedicated modules. (In fact this example was chosen such that the network creation is very simple.)
* Line 11 is just a quick way to generate the network described above. For more complex networks, you will either have to write more complex function or use dedicated modules. (In fact this example was chosen such that the network creation is very simple.)

* The values of :math:`τ` are initialised globally (line 10). We should not just define a function here, because 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 disastrous).
* The values of :math:`τ` are initialised globally (line 12). We should not just define a function here, because 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 disastrous).

* In line 15, we use `symengine.sin` – in contrast to `math.sin` or `numpy.sin`.
* In line 17, we use `symengine.sin` – in contrast to `math.sin` or `numpy.sin`.

* In line 20, we explicitly specify the delays to speed things up a little.
* In line 22, we explicitly specify the delays to speed things up a little.

* In line 21, we explicitly use absolute instead of relative errors, as the latter make no sense for Kuramoto oscillators.
* In line 23, we explicitly use absolute instead of relative errors, as the latter make no sense for Kuramoto oscillators.

* In line 24, we integrate blindly with a maximum time step of 0.1 up to the maximal delay to ensure that initial discontinuities have smoothened out.
* In line 26, we integrate blindly with a maximum time step of 0.1 up to the maximal delay to ensure that initial discontinuities have smoothened out.
"""

if __name__ == "__main__":
Expand All @@ -39,7 +39,6 @@

from jitcdde import jitcdde, t, y


rng = np.random.default_rng(seed=42)
n = 100
ω = 1
Expand Down
19 changes: 10 additions & 9 deletions examples/mackey_glass.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@

.. literalinclude:: ../examples/mackey_glass.py
:dedent: 1
:lines: 60-66
:lines: 62-67

Amongst our imports were the symbols for the state (`y`) and time (`t`), which have to be used to write down the differential equation such that JiTCDDE can process it.
Using them, we can write down the right-hand side of the differential equation as a list of expressions.
Expand All @@ -21,42 +21,43 @@

.. literalinclude:: ../examples/mackey_glass.py
:dedent: 1
:lines: 68-69
:lines: 69-70

We want the initial condition and past to be :math:`y(t<0) = 1`.
Hence we can use `constant_past`.
This automatically results in the integration starting at :math:`t=0`.

.. literalinclude:: ../examples/mackey_glass.py
:dedent: 1
:lines: 71
:lines: 72

If we calculate the derivative from our initial conditions, we obtain :math:`f(t=0) = 0.025`, which does not agree with the :math:`\\dot{y}(t=0) = 0` as explicitly defined in the initial conditions. Practically, this would result in an error if we started integrating without further precautions.
`step_on_discontinuities` makes some tailored integration steps to avoid this problem and to allow for the discontinuity to be smoothed out by temporal evolution.
(See `discontinuities` for alternatives and more details on this).

.. literalinclude:: ../examples/mackey_glass.py
:dedent: 1
:lines: 73
:lines: 74

Finally, we can perform the actual integration.
In our case, we integrate for 10000 time units with a sampling rate of 10 time units. We query the current time of the integrator (`DDE.t`) to start wherever `step_on_discontinuities` ended. `integrate` returns the state after integration, which we collect in the list `data`.
Finally, we use `numpy.savetxt` to store this to the file `timeseries.dat`.

.. literalinclude:: ../examples/mackey_glass.py
:dedent: 1
:lines: 75-78
:lines: 76-79

Taking everything together, our code is:

.. literalinclude:: ../examples/mackey_glass.py
:linenos:
:dedent: 1
:lines: 60-78
:lines: 62-79
"""


if __name__ == "__main__":
import numpy
import numpy as np

from jitcdde import jitcdde, t, y

Expand All @@ -73,6 +74,6 @@
DDE.step_on_discontinuities()

data = []
for time in numpy.arange(DDE.t, DDE.t+10000, 10):
for time in np.arange(DDE.t, DDE.t+10000, 10):
data.append( DDE.integrate(time) )
numpy.savetxt("timeseries.dat", data)
np.savetxt("timeseries.dat", data)
16 changes: 8 additions & 8 deletions examples/mackey_glass_lyap.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,15 @@

.. literalinclude:: ../examples/mackey_glass_lyap.py
:dedent: 1
:lines: 17-49
:emphasize-lines: 11-12, 19-20, 22, 28, 30-33
:lines: 19-49
:emphasize-lines: 9-10, 17-18, 20, 26, 28-31
:linenos:

Note that `integrate` does not only return local Lyapunov exponents but also the length of the time interval to which they apply (which differs from the time spanned by the `integrate` command and may even be zero). This length should be used to weigh the local Lyapunov exponents for statistical processing, like in line 31.
Note that `integrate` does not only return local Lyapunov exponents but also the length of the time interval to which they apply (which differs from the time spanned by the `integrate` command and may even be zero). This length should be used to weigh the local Lyapunov exponents for statistical processing, like in line 29.
"""

if __name__ == "__main__":
import numpy
import numpy as np
from scipy.stats import sem

from jitcdde import jitcdde_lyap, t, y
Expand All @@ -34,16 +34,16 @@
data = []
lyaps = []
weights = []
for time in numpy.arange(DDE.t, DDE.t+10000, 10):
for time in np.arange(DDE.t, DDE.t+10000, 10):
state, lyap, weight = DDE.integrate(time)
data.append(state)
lyaps.append(lyap)
weights.append(weight)

numpy.savetxt("timeseries.dat", data)
lyaps = numpy.vstack(lyaps)
np.savetxt("timeseries.dat", data)
lyaps = np.vstack(lyaps)

for i in range(n_lyap):
Lyap = numpy.average(lyaps[400:,i], weights=weights[400:])
Lyap = np.average(lyaps[400:,i], weights=weights[400:])
stderr = sem(lyaps[400:,i]) # Note that this only an estimate
print(f"{i+1}. Lyapunov exponent: {Lyap:.4f} +/- {stderr:.4f}")
Loading