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

Question

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[0]=-2*r
ca[-1]=-2*r

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

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

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[0]=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()

Answer 1

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[0]=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)