## PDE – Reaction Diffusion Equations

Here is a link to a SageMathCell implementation for the discrete logistic model reaction-diffusion system. I also give the code here:

```import scipy;
from scipy import integrate
from sage.plot.colors import rainbow

pts = 20 # number of spatial points
steps = 200 # time steps
deltat = 0.1 # size of time step
t = srange(0, steps*deltat, deltat) # list of time values

# Neumann matrix
NeuM = -2*identity_matrix(pts)
for k in [0..pts-2]:
NeuM[k+1,k]=1
NeuM[k, k+1]=1
NeuM[0,0]=-1
NeuM[pts-1, pts-1] = -1

# function to square a list
def square(x):
for k in [0..len(x)-1]:
x[k]=x[k]*x[k]
return(x)

# right hand side of our system
mu = 0.5 # diffusion constant
def F(x, t=0):
return(mu*NeuM*vector(x) + x - square(x))

# initial condition all zero
u0 = 
for k in [1..pts-1]:
u0 = u0 + 

# except for a small spike
u0=.01

# evolve using odeint (switch to RK4 at some point?)
u = integrate.odeint(F, u0, t)

#plot each location in corresponding color
colorz = rainbow(pts)
mainplot = line(zip(t,u.transpose()), color = colorz)
for k in[1..pts-1]:
mainplot += line(zip(t,u.transpose()[k]), color = colorz[k])
mainplot += line(zip(t,u.transpose()), color = "black")

# plot spatial profile at each time
colorz = rainbow(steps)
profile = list_plot(u, plotjoined=true, color = colorz)
for k in [1..steps-1]:
profile +=list_plot(u[k], plotjoined=true, color = colorz[k])
profile +=list_plot(u, plotjoined=true, color = "black")

# display plots for each location
mainplot.show(figsize=[4,3], ymax = 1.2)
# output profile time series as pdf
profile.save(filename="profile.pdf")
```
This entry was posted in Uncategorized. Bookmark the permalink.