import numpy as np
import matplotlib.pyplot as plt

from matplotlib import animation, rc
from IPython.display import HTML

%matplotlib inline 

# Laplacian operator
# dx is grid spacing, a is a 2d array
# roll shifts array by 1, in i or j directions 
# return Laplacian as an array
def laplacian2D(a, dx):
    return (- 4*a + np.roll(a,1,axis=0) + np.roll(a,-1,axis=0)\
            + np.roll(a,+1,axis=1) + np.roll(a,-1,axis=1)) / (dx**2)
# it is faster to use a whole array routine than to loop over all indices (in the notebook) 

# Reaction parts of differential equations for 
# three different types of different Reaction diffusion equations
# u,v are 2D arrays, alpha, beta could be numbers or arrays
# rdtype is a string specifying the model
def Reaction(u,v,alpha,beta,rdtype):
    if (rdtype == 'FN'):  # FitzHugh Nagumo model
        Ru,Rv=R_FN(u,v,alpha,beta)
        return Ru,Rv
    if (rdtype == 'Br'):  # Brusellator model
        Ru,Rv= R_Br(u,v,alpha,beta)
        return Ru,Rv
    if (rdtype == 'GS'):  # Gray-Scott model
        Ru,Rv = R_GS(u,v,alpha,beta)
        return Ru,Rv
    print('invalid reaction type')

#The FitzHugh-Nagumo Reaction Diffusion Model
# examples: Du, Dv, alpha, beta = 1, 100, 0.5, 20, dt=0.0002
def R_FN(u,v,alpha,beta): 
    Ru =  u - u**3 - v + alpha 
    Rv = (u - v)*beta
    return Ru,Rv
    
# Brusselator Reaction Diffusion model
# example: Du, Dv, alpha, beta = 2, 22, 5, 9 # dt=0.005
def R_Br(u,v,alpha,beta): 
    Ru = alpha - (beta+1)*u + (u**2)*v 
    Rv =  beta*u - (u**2)*v
    return Ru,Rv

# solution to problem #3
# The Gray Scott Reaction Diffusion model
# alpha is feed rate for u, beta is kill or drain rate
# alpha kills v
def R_GS(u,v,alpha,beta): 
    Ru =  -u*(v**2) + alpha*(1-u) 
    Rv =  u*(v**2) -(alpha+ beta)*v
    return Ru, Rv

# A different implimentation of the Laplacian operator 
# solution to problem 4
def laplacian2D_conv(a, dx):
    z = -a 
    # be careful with roll() axis=0 is probably up/down
    # axis=1 is probably left/right
    # + 1 is probably shifting in the opposite direction than you expect
    # nothing here depends on these issues because of symmetry
    roll_right = np.roll(a,+1,axis=0)  
    roll_left  = np.roll(a,-1,axis=0)
    roll_up    = np.roll(a,+1,axis=1)
    roll_down  = np.roll(a,-1,axis=1)
    roll_right_up    = np.roll(roll_right,+1,axis=1)
    roll_right_down  = np.roll(roll_right,-1,axis=1)
    roll_left_up     = np.roll(roll_left, +1,axis=1)
    roll_left_down   = np.roll(roll_left, -1,axis=1)
    z += 0.2*(roll_right + roll_left + roll_up + roll_down)
    z += 0.05*(roll_right_up + roll_right_down + roll_left_up + roll_left_down)
    return 4*z/(dx**2)
# without the factor of 4 here this would be 4 times weaker than the other Laplacian operator
# however sometimes people refer to the operation without the 4 as a Laplacian operator

# a class for initializing and computing and showing solutions to 
# a 2D reaction diffusion equation
# in this version I show both u,v images
class Update_RD(object):
    # do this when creating a class object, axis info is passed, 
    # along with parameters for the reaction diffusion equation
    def __init__(self, fig, ax, alpha=5.0,beta=9.0,Du=2.,Dv=22.,\
                 n=100,dt=0.005,dx=1,rdtype='Br',nsteps=1,\
                 Lap_fun=laplacian2D):
        self.success = 0
        self.alpha = alpha  # parameter for reaction
        self.beta = beta  # parameter for reaction
        self.Du = Du  # diffusion coefficient for u
        self.Dv = Dv  # diffusion coefficient for v
        self.n = n    # dimension of u,v arrays
        self.rdtype = rdtype  # type of reaction diffusion model, is a string
        
        self.dx = dx   # horizontal grid spacing
        self.dt = dt   # timestep
        self.nsteps = nsteps  # how many timesteps to take per display
        self.Lap_fun = Lap_fun   # which Laplacian operator to use
        
        # initial conditions for u,v arrays
        # I find that the Gray-Scott model needs these types of initial conditions
        if (rdtype=='GS'):
            self.u = np.zeros((n,n)) + 1.0
            self.v = np.zeros((n,n)) 
            for k in range(0,int(n/10)):
                i = int(n*np.random.rand())  # randomly chose locations to set
                j = int(n*np.random.rand())
                self.v[i,j]=np.random.rand()  # set to number in [0,1)
        else:  # whereas Brusellator and FitzhughNag models seem to work with random, 
            #  near zero initial conditions
            self.u = 0.05*np.random.rand(n, n)  # initial u data uniform distribution
            self.v = 0.05*np.random.rand(n, n)  # initial v data
        
        # axis info
        self.line0, = ax[0].plot([], [], 'k-')
        self.line1, = ax[1].plot([], [], 'k-')
        self.fig = fig  # passing this just in case we want it
        self.ax = ax
        self.ax[0].set_aspect('equal')  # so square plots
        self.ax[1].set_aspect('equal')
        
        # show the u,v images initial conditions
        im0=self.ax[0].imshow(self.u,cmap='jet')
        im1=self.ax[1].imshow(self.v,cmap='hsv')
        #self.fig.colorbar(im0, ax=ax[0],fraction=0.046, pad=0.04)
        #self.fig.colorbar(im1, ax=ax[1],fraction=0.046, pad=0.04)
        # I wish I could get the colorbar to work!
            
        fac_u=  self.Du*self.dt/self.dx**2  # relevant for checking cfl condition
        fac_v=  self.Dv*self.dt/self.dx**2  
        if (fac_u>1):
            print('cfl warning Du', fac_u)
        if (fac_v>1):
            print('cfl warning Dv', fac_v)
        
    # initialization for animation
    def init(self):
        lines = []
        self.success = 0
        self.line0.set_data([], [])
        self.line1.set_data([], [])
        return self.line0,self.line1,
    
    def update_nsteps(self):
    # do nsteps,  update u,v arrays
        for ii in range(self.nsteps):
            # compute diffusive terms, with Laplacian of your choice!
            du = self.Du*self.Lap_fun(self.u, self.dx)  
            dv = self.Dv*self.Lap_fun(self.v, self.dx)
            # compute reation terms
            Ru,Rv = Reaction(self.u,self.v,self.alpha,self.beta,self.rdtype)        
            # update u,v arrays
            self.u = self.u + (du + Ru)*self.dt
            self.v = self.v + (dv + Rv)*self.dt
            
    def return_uv(self):
        return self.u,self.v

    # called each display step!
    def __call__(self, i):
        # This way the plot can continuously run ?
        if i == 0:
            return self.init()

        self.success += 1  #otherwise
       
        self.update_nsteps()
        # do nsteps x timestep for each display
        
        im0=self.ax[0].imshow(self.u,cmap='jet')
        im1=self.ax[1].imshow(self.v,cmap='hsv')
        #self.fig.colorbar(im0, ax=ax[0])
        #self.fig.colorbar(im1, ax=ax[1])
        # it would be nice to have colorbars, but they screw up in the animation
        
        return self.line0, self.line1,

# set up the Brusselator model with some nice parameters
fig, ax = plt.subplots(1,2,figsize=(8,5))
ud_Br = Update_RD(fig,ax, alpha=5.0,beta=9,Du=2,Dv=22,\
                 n=100,dt=0.005,dx=1,rdtype='Br',nsteps=200)
            
anim_Br = animation.FuncAnimation(fig, ud_Br, frames=np.arange(50), init_func=ud_Br.init,
                     interval=100, blit=True)

# with matplotlib notebook this seems to run here, with matplotlib inline use HTML

# show the animation, this takes a few minutes to do the calculation
#HTML(anim_Br.to_html5_video())
HTML(anim_Br.to_jshtml())
# image on the left is u, that on the right is v

Once Loop Reflect

# problem #1 solution
# increase Delta t for Brusselator model with some parameters
fig, ax = plt.subplots(1,2,figsize=(8,5))
ud_Br2 = Update_RD(fig,ax, alpha=5.0,beta=9,Du=2,Dv=22,\
                 n=100,dt=0.010,dx=1,rdtype='Br',nsteps=30)
            
anim_Br2 = animation.FuncAnimation(fig, ud_Br2, frames=np.arange(25), init_func=ud_Br2.init,
                     interval=100, blit=True)

HTML(anim_Br2.to_html5_video())

Problem 1 answer:

Things start to screw up at about $\Delta t \gtrsim 0.01$. For the Brusselator model with parameters given above.

The larger diffusion coefficient is $D_v = 22$ and we have $\Delta x=1$. The CFL-like condition would be $$ \Delta t < \frac{ \Delta x^2}{{\rm max } (D_u,D_v)} = \frac{1}{22} \approx 0.05 $$ To order of magnitude our value of $\Delta t= 0.01$ is near this.

Problem 2 answer:

Looking at the Laplacian operator it depends on $\Delta x^{-2}$. The ony part of the differential equation that depends on $\Delta x$ is the diffusion term. This means you can rescale $\Delta x $ and at the same time $D_u, D_v$ and not change the system. If you increase $\Delta x$ you need to increase both $D_u,D_v$ by the square of the factor incresing $\Delta x$ also to achieve the same behavior. You can also vary $D_u,D_v$ to make your grid seem larger.

Making $\Delta x$ smaller does affect your choice for timestep.
The smaller the grid spacing, the smaller the needed timestep. You need a smaller timestep when the diffusion coefficients are larger.

Problem 3 answer

See the routine laplacian2D_conv() I ran the Brusselator with it, giving slightly smoother patterns.

# run the Brusselator model with the other Laplacian operator
fig, ax = plt.subplots(1,2,figsize=(8,5))
ud_Br3 = Update_RD(fig,ax, alpha=5.0,beta=9,Du=2,Dv=22,\
                 n=100,dt=0.005,dx=1,rdtype='Br',nsteps=200,\
                 Lap_fun=laplacian2D_conv)
            
anim_Br3 = animation.FuncAnimation(fig, ud_Br3, frames=np.arange(50), init_func=ud_Br3.init,
                     interval=100, blit=True)

HTML(anim_Br3.to_jshtml())

Once Loop Reflect

# problem 5 solution
# a class for initializing and computing 2D reaction diffusion equation
# alpha, beta now cover a range, they are now arrays 
class Update_grad(object):
    # do this when creating a class object, axis info is passed, 
    # along with parameters for the reaction diffusion equation
    def __init__(self, fig,ax, alpha=5.0,alpha_amp=0.1,beta=9,beta_amp=0.1,Du=2,Dv=22,\
                 n=100,dt=0.01,dx=1,rdtype='Br',nsteps=1):
        self.success = 0
        self.Du = Du  # diffusion coefficient for u
        self.Dv = Dv  # diffusion coefficient for v
        self.n = int(n)    # dimension of u,v arrays
        self.rdtype = rdtype  # type of reaction diffusion model, is a string
        
        self.dx = dx   # horizontal grid spacing
        self.dt = dt   # timestep
        self.nsteps = nsteps  # how many timesteps to take per display
        
        # initial conditions for Brussellator model  
        #  near zero but random initial conditions
        if (rdtype=='Br'):
            self.u = 0.05*np.random.rand(self.n, self.n)  # initial u data
            self.v = 0.05*np.random.rand(self.n, self.n)  # initial v data
         
        if (rdtype=='FN'):
            self.u = np.zeros((self.n,self.n)) 
            self.v = np.zeros((self.n,self.n)) 
            for k in range(0,int(n)):
                i = int(n*np.random.rand())  # randomly chose index
                j = int(n*np.random.rand())
                self.v[i,j]=np.random.rand()  # set to a random number in [0,1)
                self.v[(i+1)%self.n,j]=np.random.rand()  # set to a random number in [0,1)
                self.v[i,(j+1)%self.n]=np.random.rand()  # set to a random number in [0,1)
       
        # make alpha, beta arrays, periodic so no abrupt changes
        x_lin = np.arange(0,2*np.pi,2*np.pi/n)
        alpha_lin = alpha + alpha_amp*np.sin(x_lin)
        beta_lin = beta + beta_amp*np.sin(x_lin) 
        alpham,betam = np.meshgrid(alpha_lin, beta_lin, sparse=False, indexing='xy')
        self.alpha  = alpham  # varying horizontal 
        self.beta  = betam  # varying vertical
        
        # axis info
        self.line0, = ax[0].plot([], [], 'k-')
        self.line1, = ax[1].plot([], [], 'k-')
        self.ax = ax
        self.ax[0].set_aspect('equal')
        self.ax[1].set_aspect('equal')
        
        # show the u,v images initial conditions
        im0=self.ax[0].imshow(self.u,cmap='jet')
        im1=self.ax[1].imshow(self.v,cmap='hsv')
            
        fac_u=  self.Du*self.dt/self.dx**2  # relevant for checking cfl condition
        fac_v=  self.Dv*self.dt/self.dx**2  
        if (fac_u>1):
            print('cfl warning Du', fac_u)
        if (fac_v>1):
            print('cfl warning Dv', fac_v)
                      
    def update_nsteps(self):
    # do nsteps,  update u,v arrays
        for ii in range(self.nsteps):
            du = self.Du*laplacian2D(self.u, self.dx)  
            dv = self.Dv*laplacian2D(self.v, self.dx)
            # compute reation terms
            Ru,Rv = Reaction(self.u,self.v,self.alpha,self.beta,self.rdtype)        
            # update u,v arrays
            self.u = self.u + (du + Ru)*self.dt
            self.v = self.v + (dv + Rv)*self.dt
            
        return self.u, self.v
      
    # initialization for animation
    def init(self):
        lines = []
        self.success = 0
        self.line0.set_data([], [])
        self.line1.set_data([], [])
        return self.line0,self.line1,

    # called each display step!
    def __call__(self, i):
        # This way the plot can continuously run ?
        if i == 0:
            return self.init()

        self.success += 1  #otherwise
       
        self.update_nsteps()
        
        im0=self.ax[0].imshow(self.u,cmap='jet')
        im1=self.ax[1].imshow(self.v,cmap='hsv')
        
        return self.line0, self.line1,

# according to checks below alpha varies horizontally and beta varies vertically

# Brusellator with varying alpha and beta
fig, ax = plt.subplots(1,2,figsize=(8,5))
#alpha=5.0,beta=9,
a_mean = 3.0;  a_amp = 2.0 # horizontal
b_mean = 6.0;  b_amp = 2.0 # vertical
ud_Br4 = Update_grad(fig,ax, alpha=a_mean, alpha_amp=a_amp,\
                beta=b_mean, beta_amp = b_amp, Du=2,Dv=22,\
                 n=100,dt=0.005,dx=1,rdtype='Br',nsteps=200)
            
anim_Br4 = animation.FuncAnimation(fig, ud_Br4, frames=np.arange(50), init_func=ud_Br.init,
                     interval=100, blit=True)

# alpha varies left to right, beta varies vertically
HTML(anim_Br4.to_html5_video())  # this takes a while
#HTML(anim_Br4.to_jshtml())

fig, ax = plt.subplots(1,2,figsize=(8,5))
#alpha=5.0,beta=9,
a_mean = 5.0;  a_amp =  1.0
b_mean = 9.0;  b_amp = 1.0
ud5 = Update_grad(fig,ax, alpha=a_mean, alpha_amp=a_amp,\
                beta=b_mean, beta_amp = b_amp, Du=2,Dv=22,\
                 n=400,dt=0.0025,dx=1,rdtype='Br',nsteps=400)
            
anim5 = animation.FuncAnimation(fig, ud5, frames=np.arange(100), \
                    init_func=ud5.init,
                     interval=100, blit=True)
HTML(anim5.to_html5_video())

#ud_FN.nsteps=2000
for i in range(5):
    uarr,varr =ud_FN.update_nsteps()
    plt.imshow(uarr)

# checking which direction is right left vs up down
alpha = 5
alpha_amp = 0.1
beta = 9
beta_amp = 0.1
n=50
x_lin = np.arange(0,2*np.pi,2*np.pi/n)
alpha_lin = alpha + alpha_amp*np.sin(x_lin)
beta_lin = beta + beta_amp*np.sin(x_lin) 
alpham,betam = np.meshgrid(alpha_lin, beta_lin, sparse=False, indexing='xy')
plt.figure()
plt.imshow(betam)  # beta is varying vertically 
plt.figure()
plt.imshow(alpham)

<matplotlib.image.AxesImage at 0x11e8634e0>

alpha_lin

array([5.        , 5.01253332, 5.02486899, 5.03681246, 5.04817537,
       5.05877853, 5.06845471, 5.07705132, 5.08443279, 5.09048271,
       5.09510565, 5.09822873, 5.09980267, 5.09980267, 5.09822873,
       5.09510565, 5.09048271, 5.08443279, 5.07705132, 5.06845471,
       5.05877853, 5.04817537, 5.03681246, 5.02486899, 5.01253332,
       5.        , 4.98746668, 4.97513101, 4.96318754, 4.95182463,
       4.94122147, 4.93154529, 4.92294868, 4.91556721, 4.90951729,
       4.90489435, 4.90177127, 4.90019733, 4.90019733, 4.90177127,
       4.90489435, 4.90951729, 4.91556721, 4.92294868, 4.93154529,
       4.94122147, 4.95182463, 4.96318754, 4.97513101, 4.98746668])

HTML(anim_N.to_html5_video())

a_mean = 0.037 
a_amp =  a_mean/2 # horizontal
b_mean = 0.06  
b_amp = b_mean/8 # vertical

DDu = 0.2
DDv = DDu/2.01

fig, ax = plt.subplots(1,2,figsize=(8,5))
ud_G2 = Update_grad(fig,ax, alpha=a_mean, alpha_amp=a_amp,\
                beta=b_mean, beta_amp = b_amp, Du=DDu,Dv=DDv,\
                 n=400,dt=1,dx=1,rdtype='GS',nsteps=400)
            
anim_G2 = animation.FuncAnimation(fig, ud_G2, frames=np.arange(200), init_func=ud.init,
                     interval=100, blit=True)
HTML(anim_G2.to_html5_video())

# try the GS model again!

a_mean = 0.037 
a_amp =  a_mean/2 # horizontal
b_mean = 0.06  
b_amp = b_mean/2 # vertical

DDu = 0.2
DDv = DDu/2.01

fig, ax = plt.subplots(1,2,figsize=(8,5))
ud_GS = Update_grad(fig,ax, alpha=a_mean, alpha_amp=a_amp,\
                beta=b_mean, beta_amp = b_amp, Du=DDu,Dv=DDv,\
                 n=300,dt=1,dx=1,rdtype='GS',nsteps=300)
            
anim_GS = animation.FuncAnimation(fig, ud_GS, frames=np.arange(200), init_func=ud_GS.init,
                     interval=100, blit=True)
#HTML(anim_G.to_html5_video())

# set up the Brusselator model with some nice parameters
fig, ax = plt.subplots(1,2,figsize=(8,5))
ud_Br4 = Update_RD(fig,ax, alpha=2.0,beta=5.4,Du=0.3,Dv=0.3/8.0,\
                 n=100,dt=0.005,dx=1,rdtype='Br',nsteps=100)
            
#anim_Br4 = animation.FuncAnimation(fig, ud_Br, frames=np.arange(50),\
#           init_func=ud_Br4.init,interval=100, blit=True)
#HTML(anim_Br4.to_html5_video())

anim_Br4 = animation.FuncAnimation(fig, ud_Br4, frames=np.arange(100),\
           init_func=ud_Br4.init,interval=100, blit=True)
HTML(anim_Br4.to_html5_video())

ud_Br4.initialize_uv();  # reset the  initial conditions
ud_Br4.nsteps=1000  # increase the number of steps between outputs
for i in range(4):  
    ud_Br4.update_nsteps();  # integrate nsteps 
    uarr = ud_Br4.u   # return the images
    plt.figure()  # new figure
    plt.imshow(uarr)  # display the u image alone