# How can I numericaly solve a convection-diffusion equation with a large diffusion term?

I want to numerically solve the advection-diffusion equation: \begin{equation} u_t(x,t) + cu_x(x,t) = v u_{xx}(x,t) \end{equation} for $x \in [0,1]$ and $t \geq 0$ subject to the boundary conditions $u(0,t) = u(1,t) = 0$ and $u(x,0) = f(x)$. To compare the accuracy of the numerical solution, I will at first derive the analytical solution.

Analytical Solution The solution satisfies \begin{equation} u(x,t) = h(x,t) \exp \left( \alpha x + \beta t \right) \end{equation} If one sets $\alpha = c /2v$ and $\beta = - c^{2} / 4v$ the function $h(x,t)$ adheres to the conventional heat equation: $h_t(x,t) = vh_{xx}(x,t)$ suject to $h(0,t) = h(1,t) = 0$ and $h(x,0) = g(x)$. The solution is particularly inexpensive to calculate if one sets $g(x) = \sin ( 2\pi x)$. Then, the function $u$ is: \begin{equation} u(x,t) = \sin (2 \pi x) \exp \left( -4 \pi^2 v t - \dfrac{c^2 t}{4v} + \dfrac{cx}{2v} \right ) \end{equation}

Numerical Solution with BCTS and Crank-Nicolson To solve the advection-diffusion equation numerically, I use the BTCS and Crank-Nicolson algorithms. Define first \begin{align*} r &= \dfrac{v \Delta t}{\Delta x^{2}} \\ R &= \dfrac{c \Delta t}{\Delta x} \end{align*} The discretized version of the advection-diffusion equation adherring to the BTCS algorithm is: \begin{equation} u_{k}^{n}(1 + 2r) + u_{k-1}^{n} ( -R/2 - r) + u_{k+1}^{n} (R/2 -r) = u_{k}^{n-1} \end{equation} The differential equation for the Crank-Nicolson algorithm can be written as: \begin{align*} &u_{k}^{n} ( 1 + r) + u_{k+1}^{n} ( R/4 - r/2) + u_{k-1}^{n} ( -R/4 - r/2) = \\ & u_{k}^{n-1} ( 1 - r) + u_{k+1}^{n} ( -R/4 + r/2) + u_{k-1}^{n} ( R/4 + r/2) \end{align*}

Simulations When I solve the advection-diffusion algorithm with $c=1, v=2$ and $\Delta t = 0.001$ and $\Delta x \approx 0.0416$ the solution looks as follows: The analytical solution is both positive and negative, the numerical solution however is an inverted parabola. With a reduced strength of the diffusion term $v=1/6$ the solution is very accurate: How can I compute the numerical solution accurately over a large range of values for $v$?

I am particularly concerned about this question, as I solved the advection-diffusion equation to understand the algorithms and want to apply them to a non-linear PDE of which I don't know the analytical solution. The computer code used for the example is:

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

class ConvectionDiffusion(object):
"""
Class to construct solutions to the convection-diffusion equation
"""

def __init__(self, DELTA_T, M, V, C):
"""
Parameters:
-----------
DELTA_T: scalar(float):
The time step Delta_t

M: scalar(float):
The number of grid-points for x

V: scalar(float):
The constant for the diffusion term

C: scalar(float):
The constant for the advection term
"""
self.DELTA_T = DELTA_T
self.M = M
self.DELTA_X = 1 / (M - 1)
self.xVec = np.linspace(0, 1., num=M)
self.C = C
self.V = V
self.R = C * DELTA_T / self.DELTA_X
self.r = V * DELTA_T / self.DELTA_X**2
self.u0 = np.sin(2 * np.pi * self.xVec) * \
np.exp(self.C * self.xVec / (2 * self.V))

def Implicit(self, stencil, RHS, Sol, l_and_u, T):
"""
Solving a PDE with implicit algorithms through (banded) matrix inversion

Parameters:
-----------
stencil: fun
Function defining for each value u^{n-1} the stencil such that
the discretized PDE can be iterated forward

RHS: fun
Constructs the RHS of A x = b.

l_and_u: tupel(int)
The number of lower and upper diagonal elements used in stencil to
use linalg.solve_banded

T: scalar(float):
The time at which the solution to the PDE is expressed.

Sol: fun
Function generating next periods u^{n} from the solution in
the interior x and the boundary values

Return:
-------
initial: array_like(float):
The solution to the PDE at time T
"""
DELTA_T, initial = self.DELTA_T, self.u0
t = 0.0
while t < T:
A = stencil(initial)
b = RHS(initial)
x = linalg.solve_banded(l_and_u, A, b)
initial = Sol(x)
t = t + DELTA_T
return initial

def stencilCN(self, initial):
"""
Stencil for the Crank-Nicolson algorithm; Stencil * u^n = b

Parameters:
----------
initial: array_like(float):
The previous solution

Returns:
--------
A: banded_matrix(float):
The matrx in baned for to pass it to scipy.linalg.banded_solve
"""
A = np.zeros((3, len(initial) - 2))
R, r = self.R, self.r
A[0, 1:] = R / 4 - r / 2
A[1, :] = 1 + r
A[2, :-1] = -R / 4 - r / 2
return A

def stencil(self, initial):
"""
Stencil for the BCTS algorithm, ; Stencil * u^n = u^{n-1}

Parameters:
----------
initial: array_like(float):
The previous solution

Returns:
--------
A: banded_matrix(float):
The matrx in baned for to pass it to scipy.linalg.banded_solve
"""
A = np.zeros((3, len(initial) - 2))
R, r = self.R, self.r
A[0, 1:] = R / 2 - r  # uper diagonal
A[1, :] = 1 + 2 * r
A[2, :-1] = -R / 2 - r  # lower diagonal
return A

def RHS(self, initial):
"""
Pepare the RHS for the BCTS algorithm.

Parameters:
-----------
initial: array_like(float):
The previous solution u^{n-1}

Returns:
--------
b: array_like(float):
Only the interior values of u^{n-1}
"""
return initial[1: -1]

def RHSCN(self, initial):
"""
Pepare the RHS for the Crank-Nicolson algorithm.

Parameters:
-----------
initial: array_like(float):
The previous solution u^{n-1}

Returns:
--------
b: array_like(float):
A weighted sum of previous u^{n-1} values
"""
R, r = self.R, self.r
present = initial[1:-1] * (1 - r)
past =  initial[:-2] * (R / 4 + r / 2)
future = initial[2:] *(-R / 4 + r / 2)
b =  present + past + future
return b

def Sol(self, x):
"""

Parameters:
x: array_like(float):
The solution in the interior

Returns:
--------
uNew: array_like(float):
The entire solution
"""
uNew = np.zeros(len(x) + 2)
uNew[1:-1] = x
return uNew

def analyticalSol(self, x, T):
"""
Analytical Solution of the advection-difussion equation

Parameters:
-----------
x   array_like(float):
The state space
T:  scalar(float):
The time at which the solution is evaluated

Returns:
v:  array_like(float):
The solution for the PDE
"""
C, V = self.C, self.V
exponential = (- 4 * np.pi**2 * V * T
- C**2 * T / (4 * V) + C * x[1:-1] / (2 * V) )
interior = np.sin(2 * np.pi * x[1:-1]) * np.exp(exponential)
v = self.Sol(interior)
return v

def ComparisonSol(self, T):
"""
Calculates the analytical solution and the its numercial approximation
according to the BCTS, C-N and Explicit algorithm at time T

Parameters:
T:  scalar(float):
The time value

Returns:
--------
sol: array_like(float):
The array with all solutions

e:  array_like(float):
The normalized second error norm of each approximation
"""
xVec, DELTA_T, DELTA_X = self.xVec, self.DELTA_T, self.DELTA_X
sol, E = np.empty((len(self.u0), 3)), np.zeros(2)

## == Anaytical == ##
sol[:, 0] = self.analyticalSol(xVec, T)

## == BTCS == ##
sol[:, 1] = self.Implicit(self.stencil, self.RHS, self.Sol, (1, 1), T)
E = np.linalg.norm(sol[:, 0] - sol[:, 1]) / np.linalg.norm(sol[:, 0])

## == CN == ##
sol[:, 2] = self.Implicit(
self.stencilCN, self.RHSCN, self.Sol, (1, 1), T)
E = np.linalg.norm(sol[:, 0] - sol[:, 2]) /np.linalg.norm(sol[:, 0])
return sol, E

if __name__ == "__main__":
Simulation = ConvectionDiffusion(DELTA_T=0.001, M=50, V=1/6, C=1.)
y, e = Simulation.ComparisonSol(0.5)
x = Simulation.xVec
fig, ax1 = plt.subplots()
ax1.plot(x, y[:, 0], 'b-', label='Anaytical Solution')
ax1.set_xlabel('x')
ax1.set_ylabel('Axis Anaytical Solution', color='b')
ax1.tick_params('y', colors='b')
ax1.legend(loc=2)

ax2 = ax1.twinx()
ax2.plot(x, y[:, 1], 'r--', label='BTCS ')
ax2.plot(x, y[:, 2], 'r:', label='C-N ')
ax2.set_ylabel('Axis Numerical Solutions', color='r')
ax2.tick_params('y', colors='r')
ax2.legend()
fig.tight_layout()
plt.show()

• Have you tried rescaling your variables? – nicoguaro Jan 25 '17 at 16:37
• @nicoguaro: Thank you for your comment! Can you explain briefly in which situations rescaling facilitates the computation of a numerical solution? – fabian Jan 26 '17 at 9:31
• In general, nondimensionalization is a good idea in general. It might improve the condition number of your system. – nicoguaro Jan 31 '17 at 13:15

The solution you believe to be inaccurate is actually by far the more accurate one; you've simply plotted it in a very deceptive way. For $\nu=2$, the exact solution is actually no bigger than about $10^{-35}$ everywhere -- it's zero for all intents and purposes. Therefore the numerical solution is correct to 10 digits -- far better than the accuracy of your solution when $\nu=1/6$.