1D_FKPP_Solver.py

% matplotlib inline
#modified from 
#http://nbviewer.ipython.org/github/diogro/ode_examples/blob/master/The%20Fisher-Kolmogorov%20equation.ipynb?create=1
from numpy import *
from scipy.integrate import odeint
from matplotlib.pyplot import plot, xlabel, ylabel, figure, axes, style
style.use('ggplot')
#import time 
from timeit import default_timer as timer
import Fmodules 
import Cmodules

# the size of the spatial domain
# his is actual size, such as "kilometres"
L = 50.
a = -L
b = L
t0 = 0.0
t1 = 20.0*2

# the number of points in the grid
#grid_size = 100
NX = 401*2
NT = 101

# the integration times
ts = linspace(t0, t1, NT)
h = ts[1]-ts[0]
# the grid
xm = linspace(a,b,NX)
dx = xm[1]-xm[0]

kc = 0.5*h/(dx*dx)

# the initial condition, consisting of a small "square" in the middle
#y0 = zeros_like(grid)
#y0[grid_size//2 - 2:grid_size//2 + 2] = 0.1
#u0 = zeros_like(xm) # population density
#u0[NX//2 - 2 : NX//2 + 2] = 0.1 #initial values
#u0[abs(xm)<1]=1 #initial values 

def AblowitzZeppetella(x,t):
    return 1./((1+exp(-5/6*t+sqrt(6)*xm/6))**2)

u0 = AblowitzZeppetella(xm,ts[0])

def fkpp(u, t, dx, d):
    # spatial second derivative
    d2x = -2 * u
    d2x[1:-1] += u[2:] + u[:-2]
    d2x[0] += u[1]
    d2x[-1] += u[-2]
    d2x = d2x/(dx*dx)
    #add the reaction terms
    du = u * (1. - u) + d2x
    return du
    
u = odeint(fkpp, u0, ts, (dx,0))

u_anal=zeros_like(u)
for i,t in enumerate(ts):
    u_anal[i,:] = AblowitzZeppetella(xm,t)
    
from JSAnimation.IPython_display import display_animation
from matplotlib import animation

fig = figure(figsize=(8,5))
ax = axes(xlim=(-L,L), ylim=(0,1.01))
line = ax.plot([], [], color='#003366', ls='--', lw=3)[0]
lineAnal = ax.plot([], [], color='red', ls='--', lw=3)[0]
# let's plot the solution
#plot(xm, u[0,:])
xlabel('space')
ylabel('u')


def fisher_kolmogorov(i):
    line.set_data(xm,u[i,:])
    lineAnal.set_data(xm,u_anal[i,:])
    
anim=animation.FuncAnimation(fig, fisher_kolmogorov,frames=NT, interval=100)
anim.save('FKPP.gif', writer='imagemagick', fps=4);