Turing pattern using FDM.

Script is based on this webpage: https://blogs.mathworks.com/graphics/2015/03/16/how-the-tiger-got-its-stripes/

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
import time

# periodic boundary condition laplacian FD scheme
def my_laplacian(in_M):
    out = -4*in_M + (np.roll(in_M,1,axis=1) + np.roll(in_M,-1,axis=1) + np.roll(in_M,1,axis=0) + np.roll(in_M,-1,axis=0))

    return out

# model parameters
da      =  1/4.0   # diffusion A - there is factor four between the matlab script and standard FD Laplace implementation
db      = .5/4.0  # diffusion B
f       = .055 # feed rate
k       = .062 # kill rate

# Mesh and initial conditions
nx      = 128
A       = np.ones((nx,nx))
B       = np.zeros((nx,nx))
B[np.ix_(np.arange(50,61),np.arange(50,71))] = 1
B[np.ix_(np.arange(60,81),np.arange(70,81))] = 1
Anew    = A.copy()
Bnew    = B.copy()
x = np.linspace(0, 127, nx)
y = np.linspace(0, 127, nx)
xv, yv = np.meshgrid(x, y)

# time stepping
dt      = .25
tottime = 4000
nt      = int(tottime/dt)
t       = 0

# hide
# solver
def turing_solve(A, B):
    Anew = A + (da*my_laplacian(A) - A*np.square(B) + f*(1-A))*dt
    Bnew = B + (db*my_laplacian(B) + A*np.square(B) - (k+f)*B)*dt

    return Anew,Bnew

# solver
#def turing_solve(A, B):
#    Anew =???
#    Bnew = ???

#    return Anew,Bnew

# based on this solution: https://stackoverflow.com/questions/42386372/increase-the-speed-of-redrawing-contour-plot-in-matplotlib/42398244#42398244
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(10,5))
ax  = plt.axes(xlim=(0, nx), ylim=(0, nx))

timeLabel=ax.text(0.02,0.98,'Time: ',transform=ax.transAxes,va='top')
ax.set_xlabel('X')
ax.set_ylabel('Y')
cmap=plt.cm.gist_yarg
p = [ax.contour(xv,yv,B, cmap=cmap ) ]
t = np.ones(nt)*time.time()

skipf = 20

# Animation function which updates figure data.  This is called sequentially
def animate(i):
    timeLabel._text='Time: %.1f'%((i*skipf)*dt)

    # use global keyword to store the latest solution
    global Anew, Bnew

    for tp in p[0].collections:
        tp.remove()
    p[0] = ax.contourf(xv,yv,Bnew, cmap= cmap)
    t[1:] = t[0:-1]
    t[0] = time.time()

    # the movie file gets too large if we save every time step, so we always skip 10
    for it in range(0,skipf):
        Anew,Bnew=turing_solve(Anew, Bnew)

    return p[0].collections+[timeLabel]


# Call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate,
                               frames=int(nt/skipf), interval=30, repeat=False, blit=True)
plt.close(anim._fig)

# Call function to display the animation
#HTML(anim.to_html5_video())  # lower resolution
HTML(anim.to_jshtml())  # higher resolution

Once Loop Reflect