Numerical Solvers for Fractional Differential Equations#

@Chaoming Wang

import brainpy as bp
import brainpy.math as bm

import matplotlib.pyplot as plt

Factional differential equations have several definitions. It can be defined in a variety of different ways that do often do not all lead to the same result even for smooth functions. In neuroscience, we usually use the following two definitions:

  • Grünwald-Letnikov derivative

  • Caputo fractional derivative

See Fractional calculus - Wikipedia for more details.

Methods for Caputo FDEs#

For a given fractional differential equation

\[ \frac{d^{\alpha} x}{d t^{\alpha}}=F(x, t) \]

where the fractional order \(0<\alpha\le 1\). BrainPy provides two kinds of methods:

  • Euler method - brainpy.fde.CaputoEuler

  • L1 schema integration - brainpy.fde.CaputoL1Schema

brainpy.fde.CaputoEuler#

brainpy.fed.CaputoEuler provides one-step Euler method for integrating Caputo fractional differential equations.

Given a fractional-order Qi chaotic system

\[\begin{split} \left\{\begin{array}{l} D^{\alpha} x_{1}=a\left(x_{1}-x_{2}\right)+x_{2} x_{3} \\ D^{\alpha} x_{2}=c x_{1}-x_{2}-x_{1} x_{3} \\ D^{\alpha} x_{3}=x_{1} x_{2}-b x_{3} \end{array}\right. \end{split}\]

we can solve the equation system by:

a, b, c = 35, 8/3, 80

def qi_system(x, y, z, t):
    dx = -a*x + a*y + y*z
    dy = c*x - y - x*z
    dz = -b*z + x*y
    return dx, dy, dz
dt = 0.001
duration = 50
inits = [0.1, 0.2, 0.3]

# The numerical integration of FDE need to know all
# history information, therefore, we need provide
# the overall simulation time "num_step" to save
# all history values.
integrator = bp.fde.CaputoEuler(qi_system,
                                alpha=0.98,  # fractional order
                                num_memory=int(duration/dt),
                                inits=inits)

runner = bp.integrators.IntegratorRunner(integrator,
                                         monitors=list('xyz'),
                                         inits=inits,
                                         dt=dt)
runner.run(duration)
WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.x, runner.mon.y)
plt.show()
../_images/f95d9de76c81a6572ee75a67e9a656604aa0da2d36a24467a7ab77b0de14bfbe.png
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.x, runner.mon.z)
plt.show()
../_images/910e2f23eb5aca7388f5db5e36f5c1ba1e1c26e9f5096802d7acdeea6b9d46b9.png

brainpy.fde.CaputoL1Schema#

brainpy.fed.CaputoL1Schema is another commonly used method to integrate Caputo derivative equations. Let’s try it with a fractional-order Lorenz system, which is given by:

\[\begin{split} \left\{\begin{array}{l} D^{\alpha} x=a\left(y-x\right) \\ D^{\alpha} y= x * (b - z) - y \\ D^{\alpha} z =x * y - c * z \end{array}\right. \end{split}\]
a, b, c = 10, 28, 8 / 3

def lorenz_system(x, y, z, t):
    dx = a * (y - x)
    dy = x * (b - z) - y
    dz = x * y - c * z
    return dx, dy, dz
dt = 0.001
duration = 50
inits = [1, 2, 3]

integrator = bp.fde.CaputoL1Schema(lorenz_system,
                                   alpha=0.99,  # fractional order
                                   num_memory=int(duration/dt),
                                   inits=inits)

runner = bp.integrators.IntegratorRunner(integrator,
                                         monitors=list('xyz'),
                                         inits=inits,
                                         dt=dt)
runner.run(duration)
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.x, runner.mon.y)
plt.show()
../_images/60f7408bffdfc5111c92580e598c05455082d516406e9952411c7b060400efb3.png
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.x, runner.mon.z)
plt.show()
../_images/bf44ff13321d76b798f263ecc1d1eefbf9983f3eb019737d5feaa378fc106bcb.png

Methods for Grünwald-Letnikov FDEs#

Grünwald-Letnikov FDE is another commonly-used type in neuroscience. Here, we provide a efficient computation method according to the short-memory principle in Grünwald-Letnikov method.

brainpy.fde.GLShortMemory#

brainpy.fde.GLShortMemory is highly efficient, because it does not require infinity memory length for numerical solution. Due to the decay property of the coefficients, brainpy.fde.GLShortMemory implements a limited memory length to reduce the computational time. Specifically, it only relies on the memory window of num_memory length. With the increasing width of memory window, the accuracy of numerical approximation will increase.

Here, we demonstrate it by using a fractional-order Chua system, which is defined as

\[\begin{split} \left\{\begin{array}{l} D^{\alpha_{1}} x=a\{y- (1+m_1) x-0.5*(m_0-m_1)*(|x+1|-|x-1|)\} \\ D^{\alpha_{2}} y=x-y+z \\ D^{\alpha_{3}} z=-b y-c z \end{array}\right. \end{split}\]
a, b, c = 10.725, 10.593, 0.268
m0, m1 = -1.1726, -0.7872

def chua_system(x, y, z, t):
    f = m1*x+0.5*(m0-m1)*(abs(x+1)-abs(x-1))
    dx = a*(y-x-f)
    dy = x - y + z
    dz = -b*y - c*z
    return dx, dy, dz
dt = 0.001
duration = 200
inits = [0.2, -0.1, 0.1]

integrator = bp.fde.GLShortMemory(chua_system,
                                  alpha=[0.93, 0.99, 0.92],
                                  num_memory=1000,
                                  inits=inits)

runner = bp.integrators.IntegratorRunner(integrator,
                                         monitors=list('xyz'),
                                         inits=inits,
                                         dt=dt)
runner.run(duration)
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.x, runner.mon.z)
plt.show()
../_images/e540ca614b30533afb874d0004f52bbd136e804cb47055353a9e82720bf922ea.png
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.y, runner.mon.z)
plt.show()
../_images/d4f35745b711c56911206cbaceb476ce92a7a115d370a2c32151423c39137d0e.png

Actually, the coefficient used in brainpy.fde.GLWithMemory can be inspected through:

plt.figure(figsize=(10, 6))
coef = integrator.binomial_coef
alphas = bm.as_numpy(integrator.alpha)

plt.subplot(211)
for i in range(3):
    plt.plot(coef[:, i], label=r'$\alpha$=' + str(alphas[i]))
plt.legend()
plt.subplot(212)
for i in range(3):
    plt.plot(coef[:10, i], label=r'$\alpha$=' + str(alphas[i]))
plt.legend()
plt.show()
../_images/fe2131e6d2e8caad99147e218dfb546631034522edc37cf8ede2bb968afabd67.png

As you see, the coefficients decay very quickly!

Further reading#

More examples of how to use numerical solvers of fractional differential equations defined in BrainPy, please see: