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.

import brainpy as bp

bp.__version__

'1.0.3'

Explicit Runge-Kutta methods

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:

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

Adaptive Runge-Kutta methods

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:

import numpy as np
import matplotlib.pyplot as plt

# 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

times = np.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 = np.array(hist_x)
hist_y = np.array(hist_y)
hist_z = np.array(hist_z)
hist_dt = np.array(hist_dt)

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()

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

# 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

lorenz_non_adaptive(x=1., y=1., z=1., t=0., sigma=10, beta=8/3, rho=28, dt=dt)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-32a7a847de20> in <module>
      8     return dx, dy, dz
      9
---> 10 lorenz_non_adaptive(x=1., y=1., z=1., t=0., sigma=10, beta=8/3, rho=28, dt=dt)

TypeError: ode_brainpy_intg_of_lorenz_non_adaptive() got an unexpected keyword argument 'dt'

Exponential Euler methods

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 \quad\quad+ \alpha _{m}(V) \\ {\frac {dh}{dt}} &= \left[-\alpha _{h}(V)-\beta _{h}(V)\right]h \quad\quad+ \alpha _{h}(V) \\ {\frac {dn}{dt}} &= \left[-\alpha _{n}(V)-\beta _{n}(V)\right]n \quad\quad+ \alpha _{n}(V) \\ \end{aligned}\end{split}\]

Iext=10.;   ENa=50.;   EK=-77.;   EL=-54.387
C=1.0;      gNa=120.;  gK=36.;    gL=0.03

def derivative(V, m, h, n, t, Iext, gNa, ENa, gK, EK, gL, EL, C):
    alpha = 0.1 * (V + 40) / (1 - bp.ops.exp(-(V + 40) / 10))
    beta = 4.0 * bp.ops.exp(-(V + 65) / 18)
    dmdt = alpha * (1 - m) - beta * m

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

    alpha = 0.01 * (V + 55) / (1 - bp.ops.exp(-(V + 55) / 10))
    beta = 0.125 * bp.ops.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

def run(method, Iext=10.):
    hist_times = np.arange(0, 100, method.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

int1 = bp.odeint(f=derivative, method='euler', dt=0.1)

run(int1, Iext=10)

<ipython-input-9-e2231649300f>:2: RuntimeWarning: overflow encountered in exp
  alpha = 0.1 * (V + 40) / (1 - bp.ops.exp(-(V + 40) / 10))
<ipython-input-9-e2231649300f>:3: RuntimeWarning: overflow encountered in exp
  beta = 4.0 * bp.ops.exp(-(V + 65) / 18)
<ipython-input-9-e2231649300f>:6: RuntimeWarning: overflow encountered in exp
  alpha = 0.07 * bp.ops.exp(-(V + 65) / 20.)
<ipython-input-9-e2231649300f>:7: RuntimeWarning: overflow encountered in exp
  beta = 1 / (1 + bp.ops.exp(-(V + 35) / 10))
<ipython-input-9-e2231649300f>:10: RuntimeWarning: overflow encountered in exp
  alpha = 0.01 * (V + 55) / (1 - bp.ops.exp(-(V + 55) / 10))
<ipython-input-9-e2231649300f>:11: RuntimeWarning: overflow encountered in exp
  beta = 0.125 * bp.ops.exp(-(V + 65) / 80)
<ipython-input-9-e2231649300f>:4: RuntimeWarning: invalid value encountered in double_scalars
  dmdt = alpha * (1 - m) - beta * m
<ipython-input-9-e2231649300f>:8: RuntimeWarning: invalid value encountered in double_scalars
  dhdt = alpha * (1 - h) - beta * h
<ipython-input-9-e2231649300f>:12: RuntimeWarning: invalid value encountered in double_scalars
  dndt = alpha * (1 - n) - beta * n

int2 = bp.odeint(f=derivative, method='euler', dt=0.02)

run(int2, Iext=10)

RK4 Method

int3 = bp.odeint(f=derivative, method='rk4', dt=0.1)

run(int3, Iext=10)

int4 = bp.odeint(f=derivative, method='rk4', dt=0.2)

run(int4, Iext=10)

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

Exponential Euler Method

int5 = bp.odeint(f=derivative, method='exponential_euler', dt=0.2)

run(int5, Iext=10)

---

Author:

• Chaoming Wang

• Email: adaduo@outlook.com

• Date: 2021.05.29

---