Inverse Laplace Transform of $\frac{1}{s} - \frac{a}{s^2 + a s^{3/2} \coth\sqrt{s}}$

Question

I got a problem for inverse Laplace transform when solving a PDE, the solution in Laplace space is

$$ \widehat{f}(s) = \frac{1}{s} - \frac{a}{s^2 + a s^{3/2}\coth{\sqrt{s}}} $$

where $a$ is a parameter and $a>0$.

$\widehat{f}(s)$ has infinite roots on the negative part of the real axis and depend on $a$.

(1) Is that possible to work out the inversion, $ f(t) = \mathscr{L}^{-1} \{ \widehat{f}(s) \}$ ?

Intuitively, I do an expansion for $\coth\sqrt{s}$ at $s = 0$,

$$ \coth \sqrt{s} = \frac{1}{\sqrt{s}} + \frac{\sqrt{s}}{3} + O(s^{3/2}), $$

which gives an approximate form of $\widehat{f}(s)$ and inversion leads to

$$ f(t) = \frac{a-1}{a} + \frac{\exp({-\frac{3at}{a+3}})}{a}. $$

This approximate form actually describes the asymptotic behavior when $t$ is large.

(2) Is that possible to give a approximate form when $t$ is small? I read a book about inverse Laplace transform, and it seems that the steepest descent method can be used for asymptotic expansion like $$ I = \int_{\mathscr{C}} h(s)e^{\tau{g(s)}} ds, $$ when $\tau \rightarrow \infty$, which might be a way, but I don't know how to do this.

Answer 1

We want to calculate the inverse Laplace transform of $$F(s)=\frac{1}{s}-\frac{a}{g(s)}$$ with $g(s)=s^2+a s\sqrt{s}\coth\sqrt{s}$. For this, we need to compute the following Bromwich integral $$f(t)=\mathcal{L}^{-1}\left\{F(s)\right\}=\frac{1}{2\pi i}\int_{S-i\infty}^{S+i\infty}F(s)e^{st}ds,$$ with $S$ such that $F(s)$ has no poles to the right of $S$. This integral can be written, with the residue theorem, as $$f(t)=\sum_n e^{s_n t} \mathrm{Res}(F,s_n),$$ where $n$ sums over all poles of $F(s)$.

First, note that $F(s)$ has no poles for $s$ real and positive.

Second, both terms of $F(s)$ have poles at $s=0$. As $\lim_{s\to0}g(s)=as+O(s^2)$, however, we find $\mathrm{Res}(F,s=0)=0$.

Last, for $s<0$, I write $i m=\sqrt{s}$ with $m$ real and positive. The poles of $F(s)$ coincide with the zeros of $g(s)$. We get \begin{align} g(s)&=0,&\Rightarrow\quad \cot m&=\frac{m}{a} \end{align} The above transcendental equation has solutions $m_n$ at $\pi n<m_n<\pi(n+1)$, see the figure. (For large $n$, the solutions lie at the left edge of the intervals.)

which I made with:

import numpy as np
import matplotlib.pyplot as plt
from mpmath import mp
from scipy import optimize

def g(m,a):
    return np.tan(m)-a/m

def Root(a,j):
    return optimize.newton(g, j*np.pi+0.5,args=(a,))

x=np.linspace(0.01,2*np.pi,100)
y=1/np.tan(x)
# Insert a NaN where the difference between successive points is positive
y[:-1][np.diff(y) > 0] = np.nan

fig, ax = plt.subplots()
ax.set_ylim(-2, 6)
ax.plot(x,y,'k',label='cot$(m)$')
ax.plot(x,x,'r',label='$m$')
ax.plot(x,x/2,'g',label='$m/2$')
ax.plot([Root(1,0)],[Root(1,0)],marker="o")
ax.plot([Root(2,0)],[Root(2,0)/2],marker="o")
ax.set(xlabel=r'$m$') 
ax.set_xticks([0., .5*np.pi, np.pi, 1.5*np.pi, 2*np.pi])
ax.set_xticklabels(labels=["$0$", r"$\frac{1}{2}\pi$",r"$\pi$", r"$\frac{3}{2}\pi$", r"$2\pi$"])
ax.legend()

To determine the residues at those poles, observe that $$\lim_{s\to s_n}g(s)=\frac{d g}{ds}\bigg|_{s=s_n}(s-s_n)$$ (where I used that $g(s_n)=0$) and thus $$\mathrm{Res}(F,s_n)=-\frac{a}{\frac{d g}{ds}\big|_{s=s_n}}$$ I find \begin{align} \frac{d g}{ds}&=2s+\frac{a}{2}\left(3\sqrt{s}\coth\sqrt{s}-s\coth^2\sqrt{s}+s\right)\\ \frac{d g}{ds}\bigg|_{s=s_n}&=\frac{1}{2 a}\left(s_n a(1+a)-s_n^2\right) \end{align} where, for the second step, I used $\coth\sqrt{s_n}=-\sqrt{s_n}/a$.

I finally find the following semi-analytic solution $$f(t)=\sum_n e^{s_n t} \frac{2a^2}{s_n^2-s_n a(1+a)}$$ with $s_n=-m_n^2$ determined numerically from $\cot m=m/a$. In the figure below you see the above semi-analytic solution with lines for $a=1$ (black) and $a=0.1$ (green). The dots represent a numerical inversion using the Talbot algorithm.

def f(t,a):
    g=0
    r=10 #sum range
    for j in range(0,r):
        sj = -Root(a,j)**2
        n1 = np.exp(sj*t)
        n2 = 2*a**2/(sj**2-sj*(a+a**2))
        g+=  n1*n2
    return g

mp.dps = 25; mp.pretty = True
def F(s,a):
    num = a
    den = s**2+a*s**(3/2)*mp.coth(mp.sqrt(s))
    return 1/s-num/den

fp01 = lambda s: F(s,0.1)
fp1 = lambda s: F(s,1)

t = np.logspace(-1,2,200)
t2 = np.logspace(-1,2,20)

fig, ax = plt.subplots()
ax.plot(t,[f(i,1)for i in t],   'k',lw=2,label='semi-analytical, $a=1$')
ax.plot(t,[f(i,0.1)for i in t], 'g',lw=2,label='$a=0.1$')
ax.plot(t2,[mp.invertlaplace(fp01,i,method='talbot') for i in t2],'g',lw=0,marker='o')
ax.plot(t2,[mp.invertlaplace(fp1,i,method='talbot') for i in t2],'k',lw=0,marker='o',label='numerical')
ax.set_xscale('log')
ax.set(xlabel=r'$t$', ylabel=r'$f(t)$') 
ax.legend()

Regarding your speculation about the late-time behavior: your $f(t)$ incorrectly approaches $1-1/a$ for $t$. Your time-scale is correct though: $\cot m=m/a$ is equivalent to $m \tan m=a$. Its first pole was approximated in Eq. (9) of https://doi.org/10.1088/2399-6528/aac0e8 through Pad'{e} approximation as $m=\sqrt{3a/(3+a)}$. With $s=-m^2$ we then find your approximate late-time relaxation timescale.

Regarding early times, it would have helped if you had posted the original PDE. I expect it to be some diffusion-type equation on a finite interval? For such problems, one can typically find early-time behavior by taking the ``boring" unperturbed boundary to infinity, see Eqs. (4)-(6) in https://arxiv.org/abs/2103.16914. Alternatively, instead of with a Laplace transformation, one could try to solve the PDE with Green's functions. The book "Heat Conduction Using Green's Functions" by Kevin D. Cole, James V. Beck, A. Haji-Sheikh, Bahman Litkouhi extensively discussed early-time solutions.