# Reaction diffusion equation solution

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.

• You should try separation of variables one more time as the PDE is linear in $M$. Nov 10, 2015 at 7:55
• Ok, I went back and did separation of variables yet again. However, when I do that, the time function ends up $e^{\lambda t}$ or $e^{-\lambda t}$, so the solution either blows up or vanishes as $t \to \infty$. This is at odds with the steady state solution when $\frac{\partial M}{\partial t}=0$, both analytically and numerically Nov 10, 2015 at 10:12
• It depends on what frame of ref do we get steady state. It may not be possible to get steady state unless a suitable frame is chosen i.e. co-ordinate ($\zeta = x-ct$). So what makes you think that we should achieve steady state other than $0$? Nov 10, 2015 at 10:19
• Well, when you set the time derivative to 0 and treat it as a BVP only varying in x, you get solutions of the form $A\sinh(\gamma(x-1))$. This agrees with the numerical solution I get. I'm starting to wonder if I'm posing the problem correctly compared to what I'm numerically solving Nov 10, 2015 at 10:41
• The problem I have is that $M'(0,t)$ should still be a function of $t$ if we have separable equation. You can't have $M(x,0) = 0$ unless you have $M(0x,0) = X(x)T(0) = 0\implies T(0) = 0$. It seems the B.Cs are confusing. Nov 10, 2015 at 11:04

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.

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}

• Very nice to see an alternate way to solve this, and a completely different application! I developed this ODE as part of my Masters thesis in mathematical biology some time ago modelling ligand-receptor interactions in a cell field - a far cry from the Seebeck effect! Jun 21, 2022 at 9:04
• can you provide a link or reference? I would like to see it :) Jun 21, 2022 at 13:16