# Growing error from a smooth initial condition for Fisher KPP equation

I'm studying the Fisker-KPP equation on the line (and in $$]0, 100[$$ numerically): $$\partial_t u = \Delta_{xx} u + u(1-u)$$ I notice a behavior I don't understand with a smooth initial condition $$u_0$$ that has the following form: u_0(x) = \left\{ \begin{aligned} &1 \quad \mbox{if} \quad |x-50| < 10 \\ &\exp{\left(\frac{1}{3^2} - \frac{1}{||x-50|-13|^2}\right)} \left(1 - \exp{\left( -\frac{1}{||x-50|-10|^2}\right)}\right) \quad \mbox{if} \quad 10 < |x-50| < 13 \\ &0 \quad \mbox{if} \quad |x-50| > 13 \end{aligned} \right. I'm using as a numerical scheme the Strang splitting, and a code taken from here:

"""
solve a scalar diffusion-reaction equation:

phi_t = kappa phi_{xx} + (1/tau) R(phi)

using operator splitting, with implicit diffusion

M. Zingale
"""

#from __future__ import print_function

import numpy as np
from scipy import linalg
from scipy.integrate import ode
#import sys
import matplotlib.pyplot as plt

def frhs(t, phi, tau):
""" reaction ODE righthand side """
return 0.25*phi*(1.0 - phi)/tau

def jac(t, phi):
return None

def react(gr, phi, tau, dt):
""" react phi through timestep dt """

phinew = gr.scratch_array()

for i in range(gr.ilo, gr.ihi+1):
with_jacobian=False)
r.set_initial_value(phi[i], 0.0).set_f_params(tau)
r.integrate(r.t+dt)
phinew[i] = r.y[0]

return phinew

def diffuse(gr, phi, kappa, dt):
""" diffuse phi implicitly (C-N) through timestep dt """

phinew = gr.scratch_array()

alpha = kappa*dt/gr.dx**2

# create the RHS of the matrix
R = phi[gr.ilo:gr.ihi+1] + \
0.5*alpha*(    phi[gr.ilo-1:gr.ihi] -
2.0*phi[gr.ilo  :gr.ihi+1] +
phi[gr.ilo+1:gr.ihi+2])

# create the diagonal, d+1 and d-1 parts of the matrix
d = (1.0 + alpha)*np.ones(gr.nx)
u = -0.5*alpha*np.ones(gr.nx)
u[0] = 0.0

l = -0.5*alpha*np.ones(gr.nx)
l[gr.nx-1] = 0.0

# set the boundary conditions by changing the matrix elements

# homogeneous neumann
d[0] = 1.0 + 0.5*alpha
d[gr.nx-1] = 1.0 + 0.5*alpha

# dirichlet
#d[0] = 1.0 + 1.5*alpha
#R[0] += alpha*0.0

#d[gr.nx-1] = 1.0 + 1.5*alpha
#R[gr.nx-1] += alpha*0.0

# solve
A = np.matrix([u,d,l])
phinew[gr.ilo:gr.ihi+1] = linalg.solve_banded((1,1), A, R)

return phinew

def est_dt(gr, kappa, tau):
""" estimate the timestep """

# use the proported flame speed
s = np.sqrt(kappa/tau)
dt = gr.dx/s
return dt

class Grid(object):

def __init__(self, nx, ng=1, xmin=0.0, xmax=1.0, vars=None):
""" grid class initialization """

self.nx = nx
self.ng = ng

self.xmin = xmin
self.xmax = xmax

self.dx = (xmax - xmin)/nx
self.x = (np.arange(nx+2*ng) + 0.5 - ng)*self.dx + xmin

self.ilo = ng
self.ihi = ng+nx-1

self.data = {}

for v in vars:
self.data[v] = np.zeros((2*ng+nx), dtype=np.float64)

def fillBC(self, var):

if not var in self.data.keys():
sys.exit("invalid variable")

vp = self.data[var]

# Neumann BCs
vp[0:self.ilo+1] = vp[self.ilo]
vp[self.ihi+1:] = vp[self.ihi]

def scratch_array(self):
return np.zeros((2*self.ng+self.nx), dtype=np.float64)

def initialize(self):
""" initial condition """

phi = self.data["phi"]
length1 = 10.
length2 = 13.
epsilon   = length2 - length1
phi[:] = np.maximum( \
np.exp( 1./epsilon**2 - 1./(np.abs(np.abs(self.x-50.)-length2))**2) * \
( 1. - np.exp( -1./(np.abs(np.abs(self.x-50.)-length1))**2)) * \
( self.x >  50.-length2 ) * \
( self.x <  50.+length2 ) \
, \
( self.x >= 50.-length1 ) * \
( self.x <= 50.+length1 ) \
)

def interpolate(x, phi, phipt):
""" find the x position corresponding to phipt """

idx = (np.where(phi >= 0.5))[0][0]
xs   = np.array([x[idx-1],   x[idx],   x[idx+1]])
phis = np.array([phi[idx-1], phi[idx], phi[idx+1]])

xpos = 0.0

for m in range(len(phis)):
# create Lagrange basis polynomial for point m
l = None
n = 0
for n in range(len(phis)):
if n == m:
continue

if l == None:
l = (phipt - phis[n])/(phis[m] - phis[n])
else:
l *= (phipt - phis[n])/(phis[m] - phis[n])

xpos += xs[m]*l

return xpos

def evolve(nx, kappa, tau, tmax, dovis=1, return_initial=0):
"""
the main evolution loop.  Evolve

phi_t = kappa phi_{xx} + (1/tau) R(phi)

from t = 0 to tmax
"""

# create the grid
gr = Grid(nx, ng=1, xmin = 0.0, xmax=100.0,
vars=["phi", "phi1", "phi2"])

# pointers to the data at various stages
phi  = gr.data["phi"]
phi1 = gr.data["phi1"]
phi2 = gr.data["phi2"]

# initialize
gr.initialize()

phi_init = phi.copy()

# runtime plotting
if dovis == 1: plt.ion()

t = 0.0
while t < tmax:

dt = est_dt(gr, kappa, tau)

if t + dt > tmax:
dt = tmax - t

# react for dt/2
phi1[:] = react(gr, phi, tau, dt/2)
gr.fillBC("phi1")

# diffuse for dt
phi2[:] = diffuse(gr, phi1, kappa, dt)
gr.fillBC("phi2")

# react for dt/2 -- this is the updated solution
phi[:] = react(gr, phi2, tau, dt/2)
gr.fillBC("phi")

t += dt

if dovis == 1:
plt.clf()
plt.plot(gr.x, phi)
plt.grid()
plt.xlim(gr.xmin,gr.xmax)
plt.ylim(0.0,1.0)
plt.title("Reaction-Diffusion, $$t = {:3.2f}$$".format(t))
plt.draw()
plt.pause(0.1)

if return_initial == 1:
return phi, gr.x, phi_init
else:
return phi, gr.x

kappa = 1.0
tau = 0.25
nx = 256
tmax1 = 1.0
phi1, x1 = evolve(nx, kappa, tau, tmax1)


As far as I can tell, the initial condition being of class $$\cal{C}^{\infty}$$, and $$1$$ being stable, the solution should remain $$1$$ where the initial condition is $$1$$. But this is not what I observe.

Is this a numerical artifact?