# Reaction diffusion numerical scheme: discrepancy between devised scheme and scheme that works

I've been writing a PDE solver for a reaction diffusion equation in one dimension. To test it I'm using the "method of manufactured solutions" (that's the term my supervisor gave it - I'm not sure if that's a common term).

The equation is:

$$\frac{\partial u}{\partial t}=D\frac{\partial^2 u}{\partial x^2}+f(x,t,u)$$

I have chosen a specific solution: $$u(x,t)=-e^{-t}x^2(2x-3L)$$ Which is the solution to the following particular equation: $$\frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}-u+De^{-t}(12x-6L)$$

I'm breaking this into two stages. I'm using forward Euler to solve the reaction part, and Crank-Nicolson to solve the diffusion part. In other words, I'm solving: $$\frac{\partial u}{\partial t}=f(x,t,u)$$ Using the numerical scheme: $$\frac{U^{k+1}_j-U^k_{j}}{\Delta t}=-U^{k}_j+De^{-t^{k}}(12x_j-6L)$$ Which gives the next time step as: $$U^{k+1}_j=U^{k}_j(1-\Delta t)+\Delta t De^{-t^{k}}(12x_j-6L)$$ (Note: superscripts represent time steps, and subscripts represent space steps, i.e. $t^{k}$ is the kth time step) A matrix equaion for this is: $$\hat{U}^{k+1}=\hat{U}^{k}(1-\Delta t)+\Delta t De^{-t^{k}}(12\hat{X}-6L)$$

And then solving: $$\frac{U^{k+1}_{j}-U^{k}_{j}}{\Delta t}=\frac{D}{2(\Delta x)^2}\left[(U^{k+1}_{j+1}-2U^{k+1}_{j}+U^{k+1}_{j-1})+(U^{k}_{j+1}-2U^{k}_{j}+U^{k}_{j-1})\right] +f(x_j,t^{k},U^{k}_{j})$$

Which boils down to this matrix equation: $$A\hat{U}^{k+1}=B\hat{U}^k+2\Delta t\hat{f}$$ (the hats represent vectors)

Apologies if this is not explicit enough, I'm trying to avoid writing out a huge matrix equation that is a pain to format. Please ask if clarification is needed.

So, up to now is what I'll call my "devised scheme". The next bit is about my "implemented scheme" that I've written up in python.

To solve this, I'm using a tridiagonal solver to solve the matrix equation: $$\hat{U}^{k+1}=A^{-1}(B\hat{U}^k+2\Delta t\hat{f})$$

Anyway, what I've done next is to replace the $\hat{f}$ with the $\hat{U}^{k+1}$ I get from the forwards Euler step. I'm not so sure if this is the right thing to do.

When I solve this matrix equation and plot the results I get: Where the dotted lines are the numerical solution and the solid lines are the exact solution.

Now, this is clearly not correct. When I change the matrix equation to this however: $$\hat{U}^{k+1}=A^{-1}(B\hat{U}^k\highlight{-}2\Delta t\hat{f})$$ and again replace the $\hat{f}$ with the $\hat{U}^{k+1}$ I get from the forwards Euler step, I get the follwing result. Much better isn't it?

So my question is, why? What am I missing, what am I doing wrong? Please just assume that the numerical method can't be changed. I can use a different way easy enough, but I need to know what is going on here to make sure I know what I'm doing in general.

EDIT:Here is my code.

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

def fwdeuler(L,dt,T,uin,D):
nx=len(uin)
unext=np.zeros(nx)
dx = L/float(nx)
for j in xrange(nx):
unext[j]= uin[j]*(1-dt)+ D*dt*np.exp(-T)*(12*j*dx-6*L)

return unext

#parameters for the simulation
xsteps=100
tsteps=200
tf=2
L=5
D=.01

dt=tf/float(tsteps)
dx=L/float(xsteps)
r=(dt*D)/(dx**2)

#the A matrix
ac=np.ones(xsteps-1)*(-r)
ca=np.ones(xsteps-1)*(-r)

#no flux BCs
ac=-2*r
ca[-1]=-2*r

b=np.ones(xsteps)*2*(1+r)

#the matrix for the tridiagonal solver
band=np.vstack((np.hstack((,ac)),np.vstack((b,np.hstack((ca,))))))

sol=np.zeros((xsteps,tsteps))

x=np.linspace(0,L,xsteps)
#initial conditions
sol[:,0]=-(x**2)*(2*x-3*L)

#creating the B matrix
a1=np.ones(xsteps)*2*(1-r)
ta2=np.ones(xsteps-1)*r
ta3=np.ones(xsteps-1)*r

ta2=2*r
ta3[-1]=2*r

poststep=np.diag(ta1)+np.diag(ta2,1)+np.diag(ta3,-1)

for i in xrange(1,tsteps):
#first step make the RHS of the equation
#bvec, as in Ax=b
bvec=np.dot(poststep,sol[:,i-1])-2*dt*fwdeuler(L,dt,i*dt,sol[:,i-1],D)

#second step solves x=A^-1*b
sol[:,i]=la.solve_banded((1,1),band,bvec)

plt.figure(num=None, figsize=(20,16), dpi=80, facecolor='w', edgecolor='k')

for i in xrange(0,tsteps,tsteps/10):
plt.plot(x,sol[:,i],'--')
plt.plot(x,-np.exp(-dt*i)*(2*x**3-3*L*x**2))
plt.show()

• Are you solving the equation in all of space i.e. $-\infty, \infty$? What are your initial conditions? – user_of_math May 11 '15 at 3:56
• I'm solving from $x=0$ to $x=L>0$. Also I should have mentioned that I've got no flux boundary conditions at these positions – Phill May 11 '15 at 3:59
• Initial conditions are $u(x,0)=x^2(2x-3L)$ – Phill May 11 '15 at 4:00
• Are you replacing $f$ with $U^{k+1}$? That doesn't sound sane. Did you mean evaluating $f$ at $U = U^{k+1}$? Are the schemes different by only the sign? – uranix May 11 '15 at 9:18
• Yeah, I'm not sure about that. I need to include the reaction part after I've solved it separately, but I don't know how. I don't have any good justification for it, but somehow this way spits out meaningful results, as you can see in the graph. And it's not because the reaction term vanishes either, if I remove it's contribution, then the results don't even approach the exact solution – Phill May 11 '15 at 10:09

So you're trying to implement following scheme $$-r U^{k+1}_{j-1} +2(1+r)U^{k+1}_j -r U^{k+1}_{j+1} = r U^{k}_{j-1} +2(1-r)U^{k}_j +r U^{k}_{j+1} + 2\Delta t f(t_n, U^k_j)$$ with $r = \frac{D\Delta t}{\Delta x^2}$. That can be implemented as $$\text{rhs} = r U^{k}_{j-1} +2(1-r)U^{k}_j +r U^{k}_{j+1} + 2\Delta t f(t_n, U^k_j)\\ \text{Solve } \quad -r U^{k+1}_{j-1} +2(1+r)U^{k+1}_j -r U^{k+1}_{j+1} = \text{rhs for } U^{k+1}.$$

That could also be written as two stage scheme (separated by reaction and diffusion processes) $$\frac{U^{k+1/2}_j - U^{k}_j}{\Delta t} = f(t_n, U_j^k)\\ \frac{U^{k+1}_j - U^{k+1/2}_j}{\Delta t} = D \Lambda \frac{U^{k}_j + U^{k+1}_j}{2}$$ where $\Lambda u_j$ is a short notation for $\frac{u_{j-1} - 2 u_j + u_{j+1}}{\Delta x^2}$. Note, if you simply add those two stages, $U^{k+1/2}$ will go away and you'll get the original $$\frac{U^{k}_j - U^{k}_j}{\Delta t} = D \Lambda \frac{U^{k}_j + U^{k+1}_j}{2} + f(t_n, U_j^k) \tag{*}$$

But your code computes somewhat different. The function fwdeuler computes $U^{k+1/2}$, but later you plug that again into right hand side, effectively solving $$\begin{gathered} \frac{U^{k+1/2}_j - U^{k}_j}{\Delta t} = f(t_n, U_j^k)\\ \frac{U^{k+1}_j - U^{k}_j}{\Delta t} = D \Lambda \frac{U^{k}_j + U^{k+1}_j}{2} + U^{k+1/2}. \end{gathered} \tag{**}$$

There's difference between $(*)$ and $(**)$. The latter just could not be right because on the left you have dimension of $[U]T^{-1}$ and on the right mix of $DL^{-2}[U] \sim [U]T^{-1}$ and just $[U]$.

So, briefly. Either rename your fwdeuler to just f and replace it with

def f(L,T,uin,D):
nx=len(uin)
unext=np.zeros(nx)
dx = L/float(nx)
for j in xrange(nx):
unext[j]= -uin[j] + D*np.exp(-T)*(12*j*dx-6*L)
...
bvec=np.dot(poststep,sol[:,i-1])+2*dt*f(L,i*dt,sol[:,i-1],D)


or implement two stage scheme like

ta1=np.ones(xsteps)*2*(-r) # Note the difference
ta2=np.ones(xsteps-1)*r
ta3=np.ones(xsteps-1)*r

ta2=2*r
ta3[-1]=2*r

poststep=np.diag(ta1)+np.diag(ta2,1)+np.diag(ta3,-1)
...
solmid = fwdeuler(L,dt,i*dt,sol[:,i-1],D) # Computing U^{k+1/2}
# rhs = r U_{j-1}^k -2r U_j^k + r U_{j+1}^k + 2 U_j^{k+1}
bvec=np.dot(poststep,sol[:,i-1]) + 2 * solmid
sol[:,i]=la.solve_banded((1,1),band,bvec)

• Fantastic, that's exactly what I was trying to do. I just couldn't couldn't work out how to do it! +1 and answer accepted. Just out of curiosity, do you have an explanation for how my method managed to match the exact solution? It was a stroke of luck that I got something to work without knowing why, but seeing that it -does- work theremust be some logcal reason! – Phill May 12 '15 at 0:39
• Actually, nevermind, I've worked out what was going on... it was a problem with how I was using ipython :/ – Phill May 12 '15 at 4:36