Reaction diffusion equation solution

Question

This has been driving me spare for the last week, and I feel pretty bad for not being able to get a solution because (at least on the face of it) it's a pretty simple equation.

I have the following reaction diffusion equation: $$\frac{\partial M}{\partial t}=d\frac{\partial^2 M}{\partial x^2}-gM$$ With: $$\frac{\partial M}{\partial x}(0,t)=-h, \quad M(1,t)=0, \quad M(x,0)=0$$ $$0 \le x \le 1, \quad t>0, \quad d,g,h \ge 0$$

I'm searching for an analytical solution. I've tried separation of variables, I've tried basic transformations, I've looked in books, but I just can't find a solution that satisfies the BCs. I've solved it numerically and (for the values I tried, relating to the larger problem I'm working on) it's a pretty boring curve, so I don't expect there to be anything crazy going on in a solution.

Please, someone put me out of my misery! I'm fully expecting that I'm missing something obvious or making a ridiculous mistake, but I'd like to see what people answer before I put my working up here.

Answer 1

The boundary condition at $x=0$ is not homogeneous. You cannot apply separation of variables directly. Let $M(x,t)=v(x,t)+\phi(x)$, where $\phi(x)$ is the steady state, that is, the unique solution of $$ d\,\phi''-g\,\phi=0,\quad\phi'(0)=-h,\quad \phi(1)=0. $$ Then $v$ satisfies the equation $$ v_t=d\,v_{xx}-g\,v $$ with boundary condition $v_x(0,t)=0$, $v(1,t)=0$ and initial value $v(x,0)=-\phi(x)$. Finally, let $v=e^{-gt}w$. Then $w$ satisfies the equation $$ w_t=d\,w_{xx} $$ with boundary condition $w_x(0,t)=0$, $w(1,t)=0$ and initial value $w(x,0)=-\phi(x)$. You can now use separation of variables to find $w$ and $$ M=e^{-gt}\,w+\phi. $$ As $t\to\infty$, $M$ approaches the steady state solution, as it should be.

Answer 2

The problem your describing was solved in Eqs. 13, 14, 15, and Appendix A of https://doi.org/10.1103/PhysRevE.99.042136, see https://arxiv.org/abs/1902.03963 for a preprint version. That paper studied the electrolyte Seebeck effect--the tendency of different types of ions to separate when subjected to a temperature difference. Your $M$ corresponds in that paper to the local ionic charge density $q$ between two flat plates held at different temperatures. I'm following that paper, where your problem is solved through Laplace inversions.

Original problem: \begin{align} \partial_t M &= d \partial_x^2 M-g M\,,\\ \partial_x M(0,t)&=-h\,,\\ M(1,t)&=0\,,\\ M(x,0)&=0\,. \end{align} Let's clean up the notation first, $M\to f h$, $t\to t'/d$, and $g\to n^2 d$, and drop primes, \begin{align} \partial_t f &= \partial_x^2 f-n^2 f\,,\\ \partial_x f(0,t)&=-1\,,\\ f(1,t)&=0\,,\\ f(x,0)&=0\,. \end{align} A Laplace transformation $\hat{f}(x,s)=\mathcal{L}\left\{f(x,t)\right\}=\int_{0}^{\infty}{\rm d}t\, e^{st}f(x,t)$ takes the above equations to \begin{align} s\hat{f}(x,s) -f(0,x)&= \partial_x^2 \hat{f}(x,s)-n^2 \hat{f}(x,s)\nonumber\\ \partial_x^2 \hat{f}(x,s) &= m^2\hat{f}(x,s)\,, \end{align} where we used that $f(0,x)=0$ and where $m^2=n^2+s$. The boundary conditions turn into \begin{align} \partial_x \hat{f}(0,s)&=-\frac{1}{s}\,,\\ \hat{f}(1,s)&=0\,. \end{align} The above ODE for $\hat{f}(x,s)$ is solved by \begin{align} \hat{f}(x,s)=-\frac{1}{sm}\frac{\sinh[m(x-1)]}{\cosh m}\,. \end{align} Determining $f(x,t)=\mathcal{L}^{-1}\left\{\hat{f}(x,s)\right\}$ requires performing an inverse Laplace transformation, which by the residue theorem amounts to \begin{align} f(x,t)&=\sum_{s\in s_{\ell}}\text{Res}\left(\hat{f}(x,s)\exp(st),s_{\ell}\right)\,, \end{align} where the poles $s_{\ell}=\{s_{0}, s_{g}, s^{\star}_{j} \}$ of $\hat{f}(x, s)$ are located at $s_{0}=0$, $s_{n}=-n^2$, and $s^{\star}_{j}=(m^{\star}_{j})^2 -n^2$ where $m^{\star}_{j}=\pm i (j-1/2)\pi\equiv\pm i\mathcal{N}_{j}$ with $j=1,2,3\ldots$ (such that $\cosh m^{\star}_{j} =0$).

The pole $s_{0}=0$ gives the steady-state solution, \begin{align} \text{Res}\left(f(x,s)\exp(s t),0\right)=-\frac{1}{n}\frac{\sinh[n(x-1)]}{\cosh n}. \end{align} For the residue of the pole at $s_{n}=-n^2$, we consider \begin{align}\label{eq:fkd} \hat{f}(x,s)\overset{s\to-n^2}=&-\frac{1}{s}\frac{(x-1)}{1}=\frac{1}{n^2}(x-1)\,. \end{align} This implies \begin{align} \text{Res}\left(\hat{f}\exp(st),s=-n^2\right)=\lim_{s\to-n^2}(s+n^2)\hat{f}(x,s)e^{-n^2t}=0. \end{align}

For the poles at $s^{\star}_{j}$ we expand \begin{align} \cosh (m)\overset{s\to s^{\star}_{j}}=& \frac{\sinh m}{2m}\bigg|_{s=s^{\star}_{j}}\left(s- s^{\star}_{j}\right)\\ \Rightarrow \frac{1}{\cosh (m)}\overset{s\to s^{\star}_{j}}=& \frac{2i(-1)^{j}m^{\star}}{s-s^{\star}_{j}}\,, \end{align} where, going to the second line we used $m(s^{\star}_{j})=\pm m^{\star}_{j}$, and $\sinh m^{\star}_{j}=i(-1)^{j+1}$. We find \begin{align} \sum_{j\ge1}\text{Res}\left(\hat{f}\exp(st),s^{\star}_{j}\right) &=-\sum_{j\ge1} \text{Res}\left(\frac{2i(-1)^{j}\sinh [m^{\star}(x-1) ] }{ s^{\star}}\frac{\exp(s t)}{s- s^{\star}},s^{\star}_{j}\right)\\ &=-2\sum_{j\ge1}\frac{(-1)^{j}\sin[\mathcal{N}_{j}(x-1)]}{n^2+\mathcal{N}_{j}^{2}}\exp{\left[-\left( n^2+\mathcal{N}_{j}^{2}\right)t\,\right]}. \end{align}

Collecting the above terms, we find \begin{align} f(x,t)&=-\frac{1}{n}\frac{\sinh[n(x-1)]}{\cosh n}-2\sum_{j\ge1}\frac{(-1)^{j}\sin[\mathcal{N}_{j}(x-1)]}{n^2+\mathcal{N}_{j}^{2}}\exp{\left[-\left( n^2+\mathcal{N}_{j}^{2}\right)t\,\right]}, \end{align} where, again $\mathcal{N}_{j}=(j-1/2)\pi$.

Let's plot this for $n=10$ and $t=5\cdot10^{-3}$ and $t=1$, and compare to a numerical Laplace inversion of $\hat{f}(x,s)$ for the same times.

I produced this with the following code:

import numpy as np
import matplotlib.pyplot as plt
from mpmath import *

def f(t,x,n):
    fs0 = - np.sinh(n*(x-1))/(n*np.cosh(n))
    fsj=0
    for j in range(1,1000):
        Nj=(j-1/2)*np.pi
        sj=-Nj**2-n**2
        fsj+=  2*(-1)**j*np.sin(Nj*(x-1))*np.exp(sj*t)/sj
    return fs0 + fsj

def F(s,x,n):
    m = mp.sqrt(n**2+s)
    return -mp.sinh(m*(x-1))/(s*m*mp.cosh(m))

mp.dps = 15
def fnum(t,x,n):
    fp = lambda s: F(s,x,n)
    return invertlaplace(fp,t,method='talbot')

x  = np.linspace(0,1,100)
x2 = np.linspace(0,1,20)

n=10
fig, ax = plt.subplots()
ax.plot(x, [f(5e-3,i,n) for i in x],    'r',lw=2,label='$t=5\cdot10^{-3}$, analytical, 1000 terms')
ax.plot(x, [f(1,i,n) for i in x],       'k',lw=2,label='$t=1$')
ax.plot(x2,[fnum(5e-3,i,n) for i in x2],'r',lw=0,marker='o',label='numerical')
ax.plot(x2,[fnum(1,i,n) for i in x2],   'k',lw=0,marker='o',)
ax.set(xlabel=r'$x$', ylabel=r'$f(x,t)$') 
ax.legend()

Returning to the original notation, \begin{align} \frac{M(x,t)}{h}&=-\frac{1}{\sqrt{g/d} }\frac{\sinh [\sqrt{g/d} (x-1)]}{\cosh \sqrt{g/d} }-2\sum_{j\ge1}\frac{(-1)^{j}\sin\!\left[\mathcal{N}_{j}(x-1)\right]}{g/d+\mathcal{N}_{j}^2}\exp{\left[-t\left(g+d\mathcal{N}_{j}^2\right]\,\right]} \end{align}