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
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
.. _numerical_integration:
Numerical integration
=====================
By default, Brian
chooses an integration method automatically, trying to solve the equations
exactly first (for linear equations) and then resorting to numerical algorithms.
It will also take care of integrating stochastic differential equations
appropriately.
Note that in some cases, the automatic choice of integration method will not be
appropriate, because of a choice of parameters that couldn't be determined in
advance. In this case, typically you will get nan (not a number) values in the
results, or large oscillations. In this case, Brian will generate a warning to
let you know, but will not raise an error.
Method choice
-------------
You will get an ``INFO`` message telling you which integration method Brian decided to use,
together with information about how much time it took to apply the integration method
to your equations. If other methods have been tried but were not applicable, you will
also see the time it took to try out those other methods. In some cases, checking
other methods (in particular the ``'exact'`` method which attempts to solve the
equations analytically) can take a considerable amount of time -- to avoid wasting
this time, you can always chose the integration method manually (see below). You
can also suppress the message by raising the log level or by explicitly suppressing
``'method_choice'`` log messages -- for details, see :doc:`../advanced/logging`.
If you prefer to chose an integration algorithm yourself, you can do so using
the ``method`` keyword for `NeuronGroup`, `Synapses`, or `SpatialNeuron`.
The complete list of available methods is the following:
* ``'exact'``: exact integration for linear equations (alternative name: ``'linear'``)
* ``'exponential_euler'``: exponential Euler integration for conditionally
linear equations
* ``'euler'``: forward Euler integration (for additive stochastic
differential equations using the Euler-Maruyama method)
* ``'rk2'``: second order Runge-Kutta method (midpoint method)
* ``'rk4'``: classical Runge-Kutta method (RK4)
* ``'heun'``: stochastic Heun method for solving Stratonovich stochastic
differential equations with non-diagonal multiplicative noise.
* ``'milstein'``: derivative-free Milstein method for solving stochastic
differential equations with diagonal multiplicative noise
.. note::
The ``'independent'`` integration method (exact integration for a system of
independent equations, where all the equations can be analytically solved
independently) should no longer be used and might be removed in future
versions of Brian.
.. note:: The following methods are still considered experimental
* ``'gsl'``: default integrator when choosing to integrate equations with
the GNU Scientific Library ODE solver: the rkf45 method. Uses an adaptable
time step by default.
* ``'gsl_rkf45'``: Runge-Kutta-Fehlberg method.
A good general-purpose integrator according to the GSL documentation. Uses an
adaptable time step by default.
* ``'gsl_rk2'``: Second order Runge-Kutta method using GSL. Uses an adaptable
time step by default.
* ``'gsl_rk4'``: Fourth order Runge-Kutta method using GSL. Uses an adaptable
time step by default.
* ``'gsl_rkck'``: Runge-Kutta Cash-Karp method using GSL. Uses an adaptable
time step by default.
* ``'gsl_rk8pd'``: Runge-Kutta Prince-Dormand method using GSL. Uses an adaptable
time step by default.
.. admonition:: The following topics are not essential for beginners.
|
Technical notes
---------------
Each class defines its own list of algorithms it tries to
apply, `NeuronGroup` and `Synapses` will use the first suitable method out of
the methods ``'exact'``, ``'euler'`` and ``'heun'`` while `SpatialNeuron`
objects will use ``'exact'``, ``'exponential_euler'``, ``'rk2'`` or ``'heun'``.
You can also define your own numerical integrators, see
:doc:`../advanced/state_update` for details.
.. _gsl_integration:
GSL stateupdaters
-----------------
The stateupdaters preceded with the gsl tag use ODE solvers defined in the GNU
Scientific Library. The benefit of using these integrators over the ones written
by Brian internally, is that they are implemented with an adaptable timestep.
Integrating with an adaptable timestep comes with two advantages:
* These methods check whether the estimated error of the solutions returned fall
within a certain error bound. For the non-gsl integrators there is currently no
such check.
* Systems no longer need to be simulated with just one time step. That is, a bigger
timestep can be chosen and the integrator will reduce the timestep when increased
accuracy is required. This is particularly useful for systems where both slow and
fast time constants coexist, as is the case with for example (networks of neurons
with) Hodgkin-Huxley equations. Note that Brian's timestep still determines the
resolution for monitors, spike timing, spike propagation etc. Hence, in a network,
the simulation error will therefore still be on the order of ``dt``. The benefit
is that short time constants occurring in equations no longer dictate the network
time step.
In addition to a choice between different integration methods, there are a few more
options that can be specified when using GSL. These options can be specified by
sending a dictionary as the ``method_options`` key upon initialization of the object
using the integrator (`NeuronGroup`, `Synapses` or `SpatialNeuron`).
The available method options are:
* ``'adaptable_timestep'``: whether or not to let GSL reduce the timestep to
achieve the accuracy defined with the ``'absolute_error'`` and
``'absolute_error_per_variable'`` options described below. If this is set to ``False``,
the timestep is determined by Brian (i.e. the ``dt`` of the respective clock is used, see :ref:`scheduling`).
If the resulted estimated error exceeds the set error bounds, the simulation
is aborted. When using cython this is reported with an `IntegrationError`.
Defaults to ``True``.
* ``'absolute_error'``: each of the methods has a way of estimating the error that
is the result of using numerical integration. You can specify the maximum size of this
error to be allowed for any of the to-be-integrated variables in base units with this
keyword. Note that giving very small values makes the simulation slow and might result
in unsuccessful integration. In the case of using the ``'absolute_error_per_variable'``
option, this is the error for variables that were not specified individually.
Defaults to 1e-6.
* ``'absolute_error_per_variable'``: specify the absolute error per variable in its own
units. Variables for which the error is not specified use the error set with the
``'absolute_error'`` option.
Defaults to None.
* ``'max_steps'``: The maximal number of steps that the integrator will take within a
single "Brian timestep" in order to reach the given error criterion. Can be set to
0 to not set any limits. Note that without limits, it can take a very long time
until the integrator figures out that it cannot reach the desired error level. This
will manifest as a simulation that appears to be stuck.
Defaults to 100.
* ``'use_last_timestep'``: with the ``'adaptable_timestep'`` option set to True, GSL tries
different time steps to find a solution that satisfies the set error bounds.
It is likely that for Brian's next time step the GSL time step
will be somewhat similar per neuron (e.g. active neurons will have a shorter GSL time step
than inactive neurons). With this option set to True, the time step GSL found to satisfy
the set error bounds is saved per neuron and given to GSL again in Brian's next time step.
This also means that the final time steps are saved in Brian's memory and can thus
be recorded with the `StateMonitor`: it can be accessed under ``'_last_timestep'``.
Note that some extra memory is required to keep track of the last time steps.
Defaults to True.
* ``'save_failed_steps'``: if ``'adaptable_timestep'`` is set to True,
each time GSL tries a time step and it results in an estimated
error that exceeds the set bounds, one is added to the ``'_failed_steps'`` variable. For
purposes of investigating what happens within GSL during an integration step, we offer
the option of saving this variable.
Defaults to False.
* ``'save_step_count'``: the same goes for the total number of GSL steps taken in a single
Brian time step: this is optionally saved in the ``'_step_count'`` variable.
Defaults to False.
Note that at the moment recording ``'_last_timestep'``, ``'_failed_steps'``, or ``'_step_count'``
requires a call to `run` (e.g. with 0 ms) to trigger the code generation process, before the
call to `StateMonitor`.
More information on the GSL ODE solver itself can be found in its
`documentation <https://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>`_.
|
