Numerical Solvers for Fractional Differential Equations#

Colab Open in Kaggle

@Chaoming Wang

import matplotlib.pyplot as plt

import brainpy as bp
import brainpy.math as bm

bp.__version__
'2.4.0'

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.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/8f8d392e8cedc71332214d982865e509444578a6947ecf9d6bca58e814e7c0f3.png
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.x, runner.mon.z)
plt.show()
../_images/41c97d19536f6f5bcc39e35a615d0fb14bf0927a7e3d15a393eb08346807a000.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.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/8856cb6a6fddd44c5c2f6d25ff1fbeb15395196093b2ffb672e1d3b600dcab00.png
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.x, runner.mon.z)
plt.show()
../_images/a6bd69e0ed489518804f59dec6e55bd0f000d921f694030cdd6b6d7aa55634f0.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.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/019b928a7d35608424fe035336620470db6c36d7b29e3d131547085e4259f10a.png
plt.figure(figsize=(10, 8))
plt.plot(runner.mon.y, runner.mon.z)
plt.show()
../_images/c511a981d04633a6439f1a15cfd800d0601812e5b08befdf6ecdc8785a4b38fa.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/05f790a23a5460192f42737d27adfa7ba8b5e1a095779094f39678947de42fad.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: