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: the code in this cell is hidden by the author

# 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

# Use imshow for fast updates and reliable blitting
im = ax.imshow(B, cmap=cmap, origin='lower', extent=(0, nx, 0, nx))

skipf = 20

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

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

    # advance the solution by several substeps to reduce frames
    for it in range(skipf):
        Anew, Bnew = turing_solve(Anew, Bnew)

    # update image
    im.set_array(Bnew)
    return [im, 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

Animation size has reached 20979886 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.

Once Loop Reflect