analytic solution for the convection-diffusion in the interval $[0,L]$

Question

I want to solve the diffusion-convection equation for a system I need to simulate in my research. Below I describe how to solve the equation for reflecting boundary conditions on both sides and a Dirac delta as an initial condition: $$\partial_t \rho(x,t) = K \partial_{xx} \rho(x,t) -V \partial_x \rho(x,t)$$ where $K>0$, $V<0$ and $x\in [0,L]$.

The boundary conditions: $$ \begin{cases} \partial_{x}\rho\left(0,t\right)-\frac{V}{K}\rho\left(0,t\right)=0\\ \partial_{x}\rho\left(L,t\right)-\frac{V}{K}\rho\left(L,t\right)=0\\ \end{cases}$$

The initial condition: $$ \rho\left(x,0\right)=\delta\left(x-x_{0}\right) $$

Here is a plot of the solution at $t=0$, with $x_0=600$ (here $L=1200$, $K=1$, $V=-\ln \phi$ for $\phi$ golden ratio):

A python notebook the generates this plot is found here.

The problems with this plot are:

1. At the time $t=0$ is doesn't look like a Dirac delta at $x_0$.
2. For times $0<t\ll \infty$ there is an expectation the solution will look like a traveling Gaussian center at $x_0 + Vt$, which is not the case.

So my question is - does the derivation I describe correct?

The diffusion-convection equation can be solved by decomposing $\rho$ as $\rho(x,t)=A(x,t)B(x,t)$ where the components satisfy: $$ \begin{cases} \partial_t A(x,t) = K \partial_{xx} A(x,t) -V \partial_x A(x,t)\\ \partial_t B(x,t) = K \partial_{xx} B(x,t) \\ \end{cases}$$

The boundary and initial conditions for $B$ are:

\begin{cases} \partial_{x}B\left(x,t\right)|_{x=0}-\frac{V}{2K}B\left(0,t\right)=0\\ \partial_{x}B\left(x,t\right)|_{x=L}-\frac{V}{2K}B\left(L,t\right)=0\\ B\left(x,0\right)=e^{-\frac{V}{2K}x}\delta\left(x-x_{0}\right) \end{cases}

The condition for $B$ to satisfy the simple diffusion is that $A(x,t)=A_0 e^{-\frac{V^2}{4K}t} e^{\frac{V}{2K}x}$. Then by separation of variables we solve for $B$:

$$ B(x,t) = C_{+}e^{\frac{V^{2}}{4K}t}e^{\frac{V}{2K}x}+\sum_{n}C_{n}e^{-K\left(\frac{n\pi}{L}\right)^{2}t}\left(\cos\left(\frac{n\pi}{L}x\right)+\frac{VL}{2n\pi K}\sin\left(\frac{n\pi}{L}x\right)\right) $$

Since $B(x,t)$ satisfies a Sturm–Liouville problem the coefficients can be determined by the multiplication of both sides with one of the eigenfunctions followed by integration over $[0,L]$:

$$ C_{n}=\frac{2 e^{-\frac{V}{2K}x_0}}{L\left(1+\left(\frac{VL}{2n\pi K}\right)^{2}\right)}\left(\cos\left(\frac{n\pi}{L}x_{0}\right)+\frac{VL}{2n\pi K}\sin\left(\frac{n\pi}{L}x_{0}\right)\right) $$

and the coefficient of the steady-state solution:

$$ C_+ = \frac{V}{K\left(e^{\frac{V}{K}L}-1\right)} $$

The combined solutions is then

\begin{equation} \rho\left(x,t\right)=\frac{V}{K\left(e^{\frac{V}{K}L}-1\right)}e^{\frac{V}{K}x}+\frac{2}{L}e^{-\frac{V^{2}}{4K}t}e^{\frac{V}{2K}\left(x-x_{0}\right)}\times\sum_{n}\frac{1}{1+\left(\frac{VL}{2n\pi K}\right)^{2}}\left(\cos\left(\frac{n\pi}{L}x_{0}\right)+\frac{VL}{2n\pi K}\sin\left(\frac{n\pi}{L}x_{0}\right)\right)e^{-K\left(\frac{n\pi}{L}\right)^{2}t}\left(\cos\left(\frac{n\pi}{L}x\right)+\frac{VL}{2n\pi K}\sin\left(\frac{n\pi}{L}x\right)\right) \end{equation}

Answer 1

@proton I believe, there is a bug in your python code. I also tried a python and found everything okay.

# MSE Question 4738335
import numpy as np
from scipy.constants import golden
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def stable_sol(x,V,L,K):
    return V/K*np.exp(V*x/K)/(np.exp(V*L/K)-1)

# Sample values
K = 1
V = -np.log(golden)
x0 = 1
L = 2
m = 200

# Compute rho as solution of convection-diffusion PDE.
# Solution by @proton
def rho(xrange, t):
    u = np.zeros(xrange.size)
    n = np.arange(1,m+1)
    i = 0
    for x in xrange:
        numer = (np.cos(n*np.pi/L*x0)+V*L/(2*n*np.pi*K)*np.sin(n*np.pi/L*x0))*np.exp(-K*(n*np.pi/L)**2*t)*(np.cos(n*np.pi/L*x)+V*L/(2*n*np.pi*K)*np.sin(n*np.pi/L*x))
        denom = 1+(V*L/(2*n*np.pi*K))**2
        u[i] = stable_sol(x,V,L,K) + 2/L*np.exp(-(V**2)/(4*K)*t)*np.exp(V/(2*K)*(x-x0))*np.sum(numer/denom)
        i += 1
        return u

# Visualize
xrange = np.linspace(0, L, 501)

fig, ax = plt.subplots(1,1)
fig.set_size_inches(5,5)

def animate(i):
    # print(i)
    t = i*0.2/100
    ax.clear()
    y = rho(xrange, t)
    ax.plot(xrange, y)
    # Set the x and y axis to display a fixed range
    ax.set_xlim([0, L])
    ax.set_ylim([0, 10])    
    ax.set_xlabel('x')
    ax.set_ylabel('rho(x,t)')
    return

ani = animation.FuncAnimation(fig, animate, interval=100, frames=100, repeat=False)
plt.show()

# ani.save(filename="./Documents/MSE_Q4738335.gif", writer="pillow")

HINT: The first frame should show the $\rho(x,0)=\delta(x-x_0)$ function. But indeed shows a $\text{si}(x-x_0)$ function. This effect occurs on filtering the $\delta(x)$ with an ideal band pass by cutting the max. frequency ('cause we can't sum until infinity).