Numerical Solvers for ODEs

BrainPy provides several numerical methods for ordinary differential equations (ODEs). Specifically, we provide explicit Runge-Kutta methods, adaptive Runge-Kutta methods, and Exponential Euler method for ODE numerical integration. For the introductionary tutorial for how to use BrainPy provided numerical solvers for ODEs, please see the document in Quickstart/Numerical Solvers.

[1]:
import brainpy as bp

bp.__version__
[1]:
'1.1.0-alpha'
[2]:
import matplotlib.pyplot as plt

%matplotlib inline

Explicit Runge-Kutta methods for ODEs

The first category of ODE numerical integration support is the explicit Runge-Kutta (RK) methods. RK methods are a huge family of numerical methods with a wide variety of trade-offs: efficiency, accuracy, stability, etc. The supported RK methods are listed in the following table:

Methods

Keywords

Euler

euler

Midpoint

midpoint

Heun’s second-order method

heun2

Ralston’s second-order method

ralston2

RK2

rk2

RK3

rk3

RK4

rk4

Heun’s third-order method

heun3

Ralston’s third-order method

ralston3

Third-order Strong Stability Preserving Runge-Kutta

ssprk3

Ralston’s fourth-order method

ralston4

Runge-Kutta 3/8-rule fourth-order method

rk4_38rule

Users can utilize these methods by specify the method option in brainpy.odeint() with their corresponding keyword. For example:

[3]:
@bp.odeint(method='rk4')
def int_v(v, t, p):
    # do something
    return v

Adaptive Runge-Kutta methods for ODEs

The second category of ODE numerical support is the adaptive RK methods. What’s different from the explicit RK methods is that adaptive methods are designed to produce an estimate of the local truncation error in a single Runge-Kutta step, then such error can be used to adaptively control the numerical step size. Specifically, if \(error > tol\), then replace \(dt\) with \(dt_{new}\) and repeat the step. Therefore, adaptive RK methods allow the varied step size. In BrainPy, the following adaptive RK methods are provided:

Methods

keywords

Runge–Kutta–Fehlberg 4(5)

rkf45

Runge–Kutta–Fehlberg 1(2)

rkf12

Dormand–Prince method

rkdp

Cash–Karp method

ck

Bogacki–Shampine method

bs

Heun–Euler method

heun_euler

In default, the above methods are not adaptive, unless users provide a keyword adaptive=True in brainpy.odeint(). When users use the adaptive RK methods for numerical integration, the instantaneously adjusted stepsize dt will be appended in the functional arguments. Moreover, the tolerance tol for stepsize adjustment can also be controlled by users. Let’s take the Lorenz system as the example:

[4]:
# adaptively adjust stepsize

@bp.odeint(method='rkf45',
           adaptive=True, # active the "adaptive" option
           tol=0.001) # set the tolerance
def lorenz(x, y, z, t, sigma, beta, rho):
    dx = sigma * (y - x)
    dy = x * (rho - z) - y
    dz = x * y - beta * z
    return dx, dy, dz
[5]:
times = bp.math.arange(0, 100, 0.01)
hist_x, hist_y, hist_z, hist_dt = [], [], [], []
x, y, z, dt = 1, 1, 1, 0.05
for t in times:
    # should provide one more argument "dt" when using the adaptive rk method
    x, y, z, dt = lorenz(x, y, z, t, sigma=10, beta=8/3, rho=28, dt=dt)
    hist_x.append(x)
    hist_y.append(y)
    hist_z.append(z)
    hist_dt.append(dt)
hist_x = bp.math.array(hist_x)
hist_y = bp.math.array(hist_y)
hist_z = bp.math.array(hist_z)
hist_dt = bp.math.array(hist_dt)
[6]:
fig = plt.figure()
ax = fig.gca(projection='3d')
plt.plot(hist_x, hist_y, hist_z)
ax.set_xlabel('x')
ax.set_xlabel('y')
ax.set_xlabel('z')

fig = plt.figure()
plt.plot(hist_dt[:100])
plt.xlabel('Step No.')
plt.ylabel('Adaptive dt')
plt.show()
<ipython-input-6-82a8e9b5ff5f>:2: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = fig.gca(projection='3d')
../_images/tutorial_intg_ode_numerical_solvers_15_1.png
../_images/tutorial_intg_ode_numerical_solvers_15_2.png

However, when adaptive=True is not set, users cannot call numerical function with the adaptively changed dt.

[7]:
# not adaptive

@bp.odeint(method='rkf45')
def lorenz_non_adaptive(x, y, z, t, sigma, beta, rho):
    dx = sigma * (y - x)
    dy = x * (rho - z) - y
    dz = x * y - beta * z
    return dx, dy, dz

try:
    lorenz_non_adaptive(x=1., y=1., z=1., t=0., sigma=10, beta=8/3, rho=28, dt=dt)
except TypeError as e:
    print("TypeError: ", e)
TypeError:  brainpy_itg_of_ode2_lorenz_non_adaptive() got an unexpected keyword argument 'dt'

Exponential Euler methods for ODEs

Finally, BrainPy provides Exponential Euler method for ODEs. For you linear ODE systems, we highly recommend you to to use Exponential Euler methods.

Methods

keywords

Exponential Euler

exponential_euler

For a linear system,

\[{dy \over dt} = A - By\]

the exponential Euler schema is given by:

\[y(t+dt) = y(t) e^{-B*dt} + {A \over B}(1 - e^{-B*dt})\]

As you can see, for such linear systems, the exponential Euler schema is nearly the exact solution.

In BrainPy, in order to automatically find out the linear part, we will utilize the SymPy to parse user defined functions. Therefore, ones need install sympy first when using exponential Euler method.

What’s interesting, the computational expensive neuron model Hodgkin–Huxley model is a linear ODE system. In the next, you will find that by using Exponential Euler method, the numerical step can be enlarged much to save the computation time.

\[\begin{split}\begin{aligned} C_{m}{\frac {d V}{dt}}&= -\left[{\bar {g}}_{\text{K}}n^{4} + {\bar {g}}_{\text{Na}}m^{3}h + {\bar {g}}_{l} \right] V +{\bar {g}}_{\text{K}}n^{4} V_{K} + {\bar {g}}_{\text{Na}}m^{3}h V_{Na} + {\bar {g}}_{l} V_{l} + I_{syn} \\ {\frac {dm}{dt}} &= \left[-\alpha _{m}(V)-\beta _{m}(V)\right]m + \alpha _{m}(V) \\ {\frac {dh}{dt}} &= \left[-\alpha _{h}(V)-\beta _{h}(V)\right]h + \alpha _{h}(V) \\ {\frac {dn}{dt}} &= \left[-\alpha _{n}(V)-\beta _{n}(V)\right]n + \alpha _{n}(V) \\ \end{aligned}\end{split}\]
[8]:
Iext=10.;   ENa=50.;   EK=-77.;   EL=-54.387
C=1.0;      gNa=120.;  gK=36.;    gL=0.03
[9]:
def derivative(V, m, h, n, t, Iext, gNa, ENa, gK, EK, gL, EL, C):
    alpha = 0.1 * (V + 40) / (1 - bp.math.exp(-(V + 40) / 10))
    beta = 4.0 * bp.math.exp(-(V + 65) / 18)
    dmdt = alpha * (1 - m) - beta * m

    alpha = 0.07 * bp.math.exp(-(V + 65) / 20.)
    beta = 1 / (1 + bp.math.exp(-(V + 35) / 10))
    dhdt = alpha * (1 - h) - beta * h

    alpha = 0.01 * (V + 55) / (1 - bp.math.exp(-(V + 55) / 10))
    beta = 0.125 * bp.math.exp(-(V + 65) / 80)
    dndt = alpha * (1 - n) - beta * n

    I_Na = (gNa * m ** 3.0 * h) * (V - ENa)
    I_K = (gK * n ** 4.0) * (V - EK)
    I_leak = gL * (V - EL)
    dVdt = (- I_Na - I_K - I_leak + Iext) / C

    return dVdt, dmdt, dhdt, dndt
[10]:
def run(method, Iext=10., dt=0.1):
    hist_times = bp.math.arange(0, 100, dt)
    hist_V, hist_m, hist_h, hist_n = [], [], [], []
    V, m, h, n = 0., 0., 0., 0.
    for t in hist_times:
        V, m, h, n = method(V, m, h, n, t, Iext, gNa, ENa, gK, EK, gL, EL, C)
        hist_V.append(V)
        hist_m.append(m)
        hist_h.append(h)
        hist_n.append(n)

    plt.subplot(211)
    plt.plot(hist_times, hist_V, label='V')
    plt.legend()
    plt.subplot(212)
    plt.plot(hist_times, hist_m, label='m')
    plt.plot(hist_times, hist_h, label='h')
    plt.plot(hist_times, hist_n, label='n')
    plt.legend()

Euler Method

[11]:
int1 = bp.odeint(f=derivative, method='euler', dt=0.1)

run(int1, Iext=10, dt=0.1)
<ipython-input-9-fb32b1d9b119>:2: RuntimeWarning: overflow encountered in exp
  alpha = 0.1 * (V + 40) / (1 - bp.math.exp(-(V + 40) / 10))
<ipython-input-9-fb32b1d9b119>:3: RuntimeWarning: overflow encountered in exp
  beta = 4.0 * bp.math.exp(-(V + 65) / 18)
<ipython-input-9-fb32b1d9b119>:6: RuntimeWarning: overflow encountered in exp
  alpha = 0.07 * bp.math.exp(-(V + 65) / 20.)
<ipython-input-9-fb32b1d9b119>:7: RuntimeWarning: overflow encountered in exp
  beta = 1 / (1 + bp.math.exp(-(V + 35) / 10))
<ipython-input-9-fb32b1d9b119>:10: RuntimeWarning: overflow encountered in exp
  alpha = 0.01 * (V + 55) / (1 - bp.math.exp(-(V + 55) / 10))
<ipython-input-9-fb32b1d9b119>:11: RuntimeWarning: overflow encountered in exp
  beta = 0.125 * bp.math.exp(-(V + 65) / 80)
<ipython-input-9-fb32b1d9b119>:4: RuntimeWarning: invalid value encountered in double_scalars
  dmdt = alpha * (1 - m) - beta * m
<ipython-input-9-fb32b1d9b119>:8: RuntimeWarning: invalid value encountered in double_scalars
  dhdt = alpha * (1 - h) - beta * h
<ipython-input-9-fb32b1d9b119>:12: RuntimeWarning: invalid value encountered in double_scalars
  dndt = alpha * (1 - n) - beta * n
../_images/tutorial_intg_ode_numerical_solvers_29_1.png
[12]:
int2 = bp.odeint(f=derivative, method='euler', dt=0.02)

run(int2, Iext=10, dt=0.02)
../_images/tutorial_intg_ode_numerical_solvers_30_0.png

RK4 Method

[13]:
int3 = bp.odeint(f=derivative, method='rk4', dt=0.1)

run(int3, Iext=10, dt=0.1)
../_images/tutorial_intg_ode_numerical_solvers_32_0.png
[14]:
int4 = bp.odeint(f=derivative, method='rk4', dt=0.2)

run(int4, Iext=10, dt=0.2)
<ipython-input-9-fb32b1d9b119>:2: RuntimeWarning: overflow encountered in exp
  alpha = 0.1 * (V + 40) / (1 - bp.math.exp(-(V + 40) / 10))
<ipython-input-9-fb32b1d9b119>:3: RuntimeWarning: overflow encountered in exp
  beta = 4.0 * bp.math.exp(-(V + 65) / 18)
<ipython-input-9-fb32b1d9b119>:6: RuntimeWarning: overflow encountered in exp
  alpha = 0.07 * bp.math.exp(-(V + 65) / 20.)
<ipython-input-9-fb32b1d9b119>:7: RuntimeWarning: overflow encountered in exp
  beta = 1 / (1 + bp.math.exp(-(V + 35) / 10))
<ipython-input-9-fb32b1d9b119>:10: RuntimeWarning: overflow encountered in exp
  alpha = 0.01 * (V + 55) / (1 - bp.math.exp(-(V + 55) / 10))
<ipython-input-9-fb32b1d9b119>:11: RuntimeWarning: overflow encountered in exp
  beta = 0.125 * bp.math.exp(-(V + 65) / 80)
<ipython-input-9-fb32b1d9b119>:4: RuntimeWarning: invalid value encountered in double_scalars
  dmdt = alpha * (1 - m) - beta * m
<ipython-input-9-fb32b1d9b119>:8: RuntimeWarning: invalid value encountered in double_scalars
  dhdt = alpha * (1 - h) - beta * h
<ipython-input-9-fb32b1d9b119>:12: RuntimeWarning: invalid value encountered in double_scalars
  dndt = alpha * (1 - n) - beta * n
<ipython-input-9-fb32b1d9b119>:17: RuntimeWarning: invalid value encountered in double_scalars
  dVdt = (- I_Na - I_K - I_leak + Iext) / C
../_images/tutorial_intg_ode_numerical_solvers_33_1.png

Exponential Euler Method

[15]:
int5 = bp.odeint(f=derivative, method='exponential_euler', dt=0.2)

run(int5, Iext=10, dt=0.2)
../_images/tutorial_intg_ode_numerical_solvers_35_0.png