Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Solving Fisher's nonlinear reaction-diffusion equation in python. Two method are used, 1) a time step method where the nonlinear reaction term is treated fully implicitly 2) a full implicit/explicit approach where a Newton iteration is used to find the solution variable at the next time step. More information http://scicomp.stackexchange.com/que…
from __future__ import division
import numpy as np
from scipy.integrate import simps
from scipy.sparse import spdiags, coo_matrix, csc_matrix, dia_matrix, dok_matrix, identity
from scipy.sparse.linalg import spsolve
from itertools import chain
import pylab
import scipy
print "scipy", scipy.__version__
print "numpy", np.__version__
# Solves the nonlinear Fisher's equation implicit/explict discretisation
# (Crank-Nicolson method) for the diffusion term and explicit for the
# reaction term using time stepping. Dirichlet BCs on l.h.s and r.h.s..
#
# u_t = d u_{xx} + B u (1 - u)
#
# BUGS: Neumann boundary conditions on the r.h.s. do not work correctly
# but provided we don't iterate the equation to forward in time
# the analytical solution can be retrieved (because the boundary
# point does not play a role in the solution).
#
# The analytical experssion is only valid for DBCs on the l.h.s.
# and NBCs on the r.h.s., but again let's not time step too far
# so that we get good agreement.
def M(x):
"""Returns coefficient matrix for diffusion equation (constant coefficients)"""
ones = np.ones(len(x))
M = spdiags( [ones, -2*ones, ones], (-1,0,1), len(x), len(x) ).tocsc()
J = len(x)
M[J-1,J-2] = 2.0 # Neumann conditions
return M
def set_DBCs(M):
M = M.todok()
M[0,0] = 1; M[0,1] = 0;
M[1,0] = 0;
M[-2,-1] = 0
M[-1,-1] = 1; M[-1,-2] = 0
return M.tocsc()
def bc_vector(x, d, h, left_value, right_value):
bc = dok_matrix((len(x), 1)) # column vector
bc[0,0] = -left_value
bc[1,0] = d/h**2 * left_value
bc[-2,0] = d/h**2 * right_value
bc[-1,0] = -right_value
return bc.tocsc()
def Imod(x):
data = np.ones(len(x))
data[0] = 0
data[-1] = 0
offset = 0
shape = (len(x),len(x))
Imod = dia_matrix((data, offset), shape=shape).todok()
return Imod.tocsc()
def solution_fishers_eqn(x,t,d,B,n):
"""Analytical solution to Fishers equation with u(x0,t)=1
and u_x(xJ,t)=0. The parameter n can be either 1 or 2."""
if n == 1:
U = 5 * np.sqrt(d * B / 6)
return 1.0 / (1 + np.exp(np.sqrt(B/(6*d))*(x - U*t)) )**2
elif n == 2:
V = np.sqrt(d*B/2)
return 1.0/ (1 + np.exp(np.sqrt(B/(2*d))*(x - V*t)) )
else:
raise ValueError("The parameter `n` must be an integer with the value 1 or 2.")
# Coefficients and simulation constants
x = np.linspace(-10, 10, 100)
x_range = np.max(x) - np.min(x)
h = x_range / len(x)
d = 0.1
B = 1.0
s = d / h**2
tau = 0.5
fsn = 1
# Initial conditions (needs to end up as a column vector)
u_init = np.matrix( solution_fishers_eqn(x,0.0,d,B,fsn) ).reshape((len(x),1))
u = u_init
# Matrices and vectors
M = M(x) # Coefficient matrix
Im = Imod(x) # Modifed identity matrix
I = identity(len(x)) # Standard identity matrix
theta = 0.5 * np.ones(len(x)) # IMEX but with a fully
theta[0] = 0.0 # implicit at boundary
theta[-1] = 0.0 # points.
theta_matrix = dia_matrix((theta, 0), shape=(len(x),len(x))).tocsc()
theta_matrix_1m = dia_matrix((1-theta, 0), shape=(len(x),len(x))).tocsc()
# Boundary conditions
#M = set_DBCs(M)
left_value = 1.0
right_value = 0.0
bc = bc_vector(x, d, h, left_value, right_value)
print "left_value", left_value
print "right_value", right_value
# Cache
slices = [x]
# Time stepping loop
STEPS = 20
PLOTS = 5
for i in range(0, STEPS+1):
u_array = np.asarray(u).flatten()
# Reaction (nonlinear, calculated everytime)
r = np.matrix(B*u_array**(fsn)*(1 - u_array)).reshape((len(x),1))
# Export data
if i % abs(STEPS/PLOTS) == 0:
plot_num = abs(i/(STEPS/PLOTS))
fraction_complete = plot_num / PLOTS
print "fraction complete", fraction_complete
print "I: %g" % (simps(u_array,dx=h), )
pylab.figure(1)
pylab.plot(x,u_array, "o", color=pylab.cm.Accent(fraction_complete))
pylab.plot(x, solution_fishers_eqn(x, i*tau, d, B, fsn), "-", color=pylab.cm.Accent(fraction_complete))
slices.append(u_array)
# Prepare l.h.s. @ t_{n+1}
A = (I - tau*s*theta_matrix*M)
# Prepare r.h.s. @ t_{n}
b = csc_matrix( (I + tau*s*theta_matrix_1m*M)*u + tau*r)
# Set boundary conditions
b[0,0] = left_value
#b[-1,-1] = right_value
# Solve linear
u = spsolve(A,b) # Returns an numpy array,
u = np.matrix(u).reshape((len(x),1)) # need a column matrix
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
np.savetxt("slices.txt", np.column_stack(slices))
pylab.title("Time-step method (dots), analytical results (lines) with dt=%g." % tau)
pylab.xlabel("x")
pylab.ylabel("u(x)")
pylab.ylim(ymax=1.2)
pylab.savefig("timestep_solution.pdf")
pylab.show()
from __future__ import division
import scipy
import numpy
print "scipy", scipy.__version__
print "numpy", numpy.__version__
import numpy as np
from scipy.sparse.linalg import inv
from scipy.integrate import simps
from scipy.sparse import spdiags, coo_matrix, csc_matrix, dia_matrix, dok_matrix, identity
from scipy.sparse.linalg import spsolve
from itertools import chain
import pylab
# Solves the nonlinear Fisher's equation implicit/explict discretisation
# (Crank-Nicolson method). A Newton iteration is used to find the unknown
# solution variable at the future time step. Dirichlet BCs on l.h.s and
# r.h.s.
#
# u_t = d u_{xx} + B u (1 - u)
#
# BUGS: Neumann boundary conditions on the r.h.s. do not work correctly
# but provided we don't iterate the equation to forward in time
# the analytical solution can be retrieved (because the boundary
# point does not play a role in the solution).
#
# The analytical experssion is only valid for DBCs on the l.h.s.
# and NBCs on the r.h.s., but again let's not time step too far
# so that we get good agreement.
def J(x,u,a,b,h,theta,n=1):
"""The Jacobian of Fishers equation with respect to the finite difference stencil."""
ones = np.ones(len(u))
u = np.asarray(u).flatten()
dFdu1j = a / h**2 * (1-theta) * np.ones(len(u))
if n == 1:
dFuj = -2*a/h**2 * (1-theta) * ones - b * u * (1 - theta) + b * (1 - theta) * (1 - u)
elif n == 2:
dFuj = -2*a/h**2 * (1-theta) - b * u * (1-theta) + b * (1-theta) * (1 - u)
else:
raise ValueError("The `n` kwargs must be an integer with the value 1 or 2.")
dFduj1 = a / h**2 * (1-theta) * np.ones(len(u))
return spdiags( [dFdu1j, dFuj, dFduj1], (-1,0,1), len(x), len(x) ).tocsc()
def M(x):
"""Returns coefficient matrix for diffusion equation (constant coefficients)
with Dirichlet BC on the left and Neumann BC on the right."""
ones = np.ones(len(x))
M = spdiags( [ones, -2*ones, ones], (-1,0,1), len(x), len(x) ).todok()
J = len(ones)
M[J-1,J-2] = 2.0 # Neumann conditions
return M.tocsc()
def solution_fishers_eqn(x,t,d,B,n):
"""Analytical solution to Fishers equation with u(x0,t)=1
and u_x(xJ,t)=0. The parameter n can be either 1 or 2."""
if n == 1:
U = 5 * np.sqrt(d * B / 6)
return 1.0 / (1 + np.exp(np.sqrt(B/(6*d))*(x - U*t)) )**2
elif n == 2:
V = np.sqrt(d*B/2)
return 1.0/ (1 + np.exp(np.sqrt(B/(2*d))*(x - V*t)) )
else:
raise ValueError("The parameter `n` must be an integer with the value 1 or 2.")
# Coefficients and simulation constants
x = np.linspace(-10, 10, 100)
x_range = np.max(x) - np.min(x)
h = x_range / len(x) # mesh spacing
d = 0.1 # diffusion coefficient
B = 1.0 # reaction constant
s = d / h**2
tau = 0.5 # time step
fsn = 1 # fishers 'n' constant (either 1 or 2)
# Boundary conditions
left_value = 1.0
right_flux = 0.0
# Initial conditions (needs to end up as a column vector)
u_init = np.matrix( solution_fishers_eqn(x,0.0,d,B,fsn) ).reshape((len(x),1))
u_init[0] = left_value # This initial value is how DBCs are incorporated the lhs is also the lhs DBC
u = u_init
# Matrices and vectors
M = M(x) # Coefficient matrix
I = identity(len(x)) # Standard identity matrix
theta = 0.5
# Cache
slices = [x]
# Time stepping loop
STEPS = 20
PLOTS = 5
for i in range(0, STEPS+1):
u_array = np.asarray(u).flatten()
if i % abs(STEPS/PLOTS) == 0:
plot_num = abs(i/(STEPS/PLOTS))
fraction_complete = plot_num / PLOTS
print "fraction complete", fraction_complete
print "I: %g" % (simps(u_array,dx=h), )
pylab.figure(1)
pylab.plot(x,u_array, "o", color=pylab.cm.Accent(fraction_complete))
pylab.plot(x, solution_fishers_eqn(x, i*tau, d, B, fsn), "-", color=pylab.cm.Accent(fraction_complete))
slices.append(u_array)
# Reaction (nonlinear, calculated everytime)
r = np.matrix(B*u_array**(fsn)*(1 - u_array)).reshape((len(x),1))
# Prepare Newton-Raphson variables
un = np.matrix(u_array).transpose()
r = np.matrix(B*np.asarray(u)**(fsn)*(1 - np.asarray(u)))
Fn = s*M*u + r
Fn[0] = 0.0 # DBC on the l.h.s.
Fn[-1] = 0.0 # DBC on r.h.s
vk = np.matrix(un + tau*Fn) # Guess value
# Modified newton iteration pg. 127 Hundsdorfer
for k in range(0,10):
r = B*np.asarray(vk)**(fsn)*(1 - np.asarray(vk))
Fk = s*M*vk + r
Fk[0] = 0.0 # DBC on the l.h.s.
Fk[-1] = 0.0 # DBC on r.h.s.
g = vk - un - (1-theta)*tau*Fn - theta*tau*Fk
J_term = I - theta*tau*J(x,un,d,B,h,theta,n=fsn)
J_term_inv = inv(J_term)
vk1 = vk - J_term_inv*g
vk = vk1
u = vk
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
np.savetxt("slices.txt", np.column_stack(slices))
pylab.title("Modified Newton solver (dots), analytical results (lines) with dt=%g." % tau)
pylab.xlabel("x")
pylab.ylabel("u(x)")
pylab.ylim(ymax=1.2)
pylab.savefig("newton_solution.pdf")
pylab.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.