# Analytical Approximations in Short Times of Exact Operational Solutions to
Reaction–Diffusion Problems on Bounded Intervals

Kwassi Anani
Department of Mathematics
University of Lomé
02 BP 1515, Maritime District, Lomé, Togo
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

This paper aims at obtaining, by means of integral transforms, analytical
approximations in short times of solutions to boundary value problems for the
one-dimensional reaction-diffusion equation with constant coefficients. The
general form of the equation is considered on a bounded generic interval and
the three classical types of boundary conditions, i.e., Dirichlet as well as
Neumann and mixed boundary conditions are considered in a unified way. The
Fourier and Laplace integral transforms are successively applied and an exact
solution is obtained in the Laplace domain. This operational solution is
proven to be the accurate Laplace transform of the infinite series obtained by
the Fourier decomposition method and presented in the literature as solutions
to this type of problem. On the basis of this unified operational solution,
four cases are distinguished where innovative formulas expressing consistent
analytical approximations in short time limits are derived with respect to the
behavior of the solution at the boundaries. Compared to the infinite series
solutions, the analytical approximations may open new perspectives and
applications, among which can be noted the improvement of numerical efficiency
in simulations of one-dimensional moving boundary problems, such as in Stefan
models.

_K_ eywords Reaction–diffusion equation; Fourier transform; Duhamel principle;
Laplace domain; Exact solutions; Asymptotic expansions; Time-step; Analytical
approximations; Stefan problems

##  1 Introduction

Initial boundary value problems for linear parabolic equations are still much
used nowadays as first model approximations of nonlinear and time-dependent
problems in bounded domains (see for examples [4] and [5]). Until now, exact
analytical solutions to such linear boundary value problems are mainly
obtained as infinite series of functions as in [20] or in [13]. But, accurate
analytical approximations in short times limits are hardly derivable from such
infinite series solutions (see [27] and [6]). The series solutions are
generally obtained via the Fourier decomposition method, by means of
separation of variables and using the regular Sturm-Liouville theory, as
presented in many textbooks of which [11] and [9] can be cited among others.
As shown in [13], any one-dimensional linear parabolic equation with constant
coefficients and a convection term can be reduced to an equation without
convection term. Similarly, the linear reaction term with a fixed coefficient
can also be canceled by specific transformations of the equation. Still, this
latter reduction is not made here for reasons of later generalization.

Compared to the classical numerical schemes as reported in [22], recent
numerical approaches to boundary value problems for heat transfer equations
have much gained in sophistication and accuracy as it can be seen in [18],
[26] and [17] for examples. However, apart of the Fourier decomposition method
mentioned above, there exist no other established procedure for finding exact
analytical solutions to boundary value problems involving linear reaction-
diffusion equations even with constant coefficients. According to [20], the
Laplace integral transform and its inversion formula are found to be not
appropriate for solving boundary value problems with a non-uniform function as
initial condition. Moreover, the classical Fourier integral applied to space
coordinates is validated only for infinite and semi-infinite solids.
Concerning other analytical and semi-approximate methods, some powerful
techniques have been used to handle heat transfer problems in finite or
infinite domains but without exhibiting new exact solutions. One can cite the
Adomian’s decomposition method as in [23] and [3], the variational iteration
method as in [19], the homotopy perturbation method as in [10] and Bessel
collocation method as in [28]. More details on analytical methods applicable
to all types of equations are available in the book by [29].

In the present paper, an exact solution in the Laplace domain with
computationally efficient analytical approximations in short time limits are
obtained to the boundary value problem for the linear reaction-diffusion
equation with constant coefficients. In section 2, the general form of the
problem is introduced and the assumptions are briefly commented. Section 3
presents the application of the Fourier integral transform and the resulting
solution is written under the form of integrals. In section 4, an explicit
solution in the Laplace domain is deduced from the previous integral form of
the solution. Then, analytical solutions in short time limits are obtained in
section 5, especially with respect to the behavior of the operational solution
at the ends points of the domain. The results are briefly applied to a
specific example in section 6 and the unified solution in the Laplace domain
has been extended to unbounded space domains. Finally, section 7 outlines the
conclusion.

##  2 The reaction–diffusion equation

We consider the one-dimensional linear parabolic equation with constant
coefficients, expressed in the nonhomogeneous reaction-diffusion form on a
bounded generic interval. For space variable $x\in[l_{1},l_{2}]$,
($l_{1},l_{2}$ reals, $l_{1}<l_{2}$) , and for time $t\in[0,T]$, ($T>0$ or
$t>0$ in the case when $T=+\infty$), the equation is stated as:

$\frac{\partial u}{\partial t}-a^{2}\frac{\partial^{2}u}{\partial
x^{2}}+bu=f(x,t),\quad l_{1}<x<l_{2},\quad 0<t<T.$ (1)

The equation (1) is subject to the initial condition

$u(x,0)=\varphi(x),\quad l_{1}\leq x\leq l_{2},$ (2)

and to the boundary conditions:

$\alpha_{1}u(l_{1},t)+\beta_{1}\frac{\partial u}{\partial
x}(l_{1},t)=g_{1}(t),\quad\alpha_{1},\beta_{1}\in\mathbb{R},\quad\alpha_{1}^{2}+\beta_{1}^{2}\neq
0,$ (3) $\alpha_{2}u(l_{2},t)+\beta_{2}\frac{\partial u}{\partial
x}(l_{2},t)=g_{2}(t),\quad\alpha_{2},\beta_{2}\in\mathbb{R},\quad\alpha_{2}^{2}+\beta_{2}^{2}\neq
0,$ (4)

for $t\in[0,T]$. The four terms in equation (1) represent respectively
transient, diffusion, reaction and source terms. The function $u(x,t)$ is to
be determined, and may represent species concentration for mass transfer or
temperature for heat transfer, while the functions $f(x,t),g_{1}(t)$ and
$g_{2}(t)$ are given. The coefficient $a>0$ is related to the constant
diffusivity of the mass or heat transfer. The reaction term (linear term in
$u$) indicates the possibility of mass or heat exchange with the environment
through the lateral surface of the body, at the rates proportional to the
concentration or to the temperature [12]. In a process of mass diffusion for
example, $b$ is the coefficient of disintegration ($b<0$) or multiplication
($b>0$). Without lost of generality, we will assume from now on that $b\geq
0$. The source term expression $f(x,t)$ may permit to consider homogeneous as
well as nonhomogeneous equations, while the function $\varphi(x)$ may indicate
zero or non-zero initial condition.

The Dirichlet, Neumann and mixed boundary conditions are expressed in a
unified way by equations (3) and (4), that is, both homogeneous and
nonhomogeneous forms of those three types of boundary conditions can be taken
into account. It is sufficient to give some acceptable values to the real
parameters $\alpha_{1},\beta_{1},\alpha_{2},\beta_{2}$, and some convenient
expressions to the time-dependent functions $g_{1}(t)$ and $g_{2}(t)$. Thus,
Dirichlet conditions are satisfied when $\alpha_{1}=\alpha_{2}=1$ and
$\beta_{1}=\beta_{2}=0$, while Neumann conditions can be obtained on the
boundaries if $\alpha_{1}=\alpha_{2}=0$ and $\beta_{1}=\beta_{2}=1$. Likewise
mixed boundary conditions are obtained when $\alpha_{1}=\alpha_{2}=1$ and
$\beta_{1}\neq 0,\ \beta_{2}\neq 0$. As it can be noted, a given combination
of these three classical types of boundary conditions can also be obtained.
Classically, the problem (1)-(4) is first solved for homogeneous equation and
boundary conditions ($f(x,t)=0$ and $g_{1}(t)=g_{2}(t)=0$), by using the
Fourier decomposition method through the separation of variables, together
with the Sturm-Liouville theory of eigenvalues and eigenfunctions. Then, the
principle of Duhamel intervenes in addition, when accounting for a non-zero
source term $f(x,t)$. Now, in the case of nonhomogeneous boundary conditions,
the problem (1)-(4) is first reduced to a problem with boundary conditions
equal to zero by the means of the so-called auxiliary functions (see [4] for
example).

The application of the above-mentioned techniques to the problem (1)-(4) will
lead to an exact series solution, which converges uniformly as well as the
series obtained by differentiating twice by $x$ and once by $t$, provided
that, $\varphi(x)$ and $f(x,t)$ are assumed continuous on $[l_{1},l_{2}]$ and
on $[l_{1},l_{2}]\times[0,T]$, and $g_{1}(t)$ and $g_{2}(t)$ are continuously
differentiable on $[0,T]$. The uniqueness of the series solution can be proven
by using the maximum principle. In general, existence and uniqueness theorems
of classical as well as weak or generalized solutions to initial-boundary-
value problems for the linear parabolic equation in one or higher dimensional
spaces, have been proven only under certain assumptions for specific classes
of functions (see [21] and [7] as examples). Let us admit here in a first step
that the source term $f(x,t)$ is continuous on $[l_{1},l_{2}]\times[0,T]$ and
the initial condition function $\varphi(x)$ satisfies at least the so-called
Dirichlet conditions relatively to the space variable $x$ on $[l_{1},l_{2}]$.
That is, this function is piece-wise continuous or can be expressed in a
unique way as a convergent series of eigenvalues and eigenfunctions of a
Sturm-Liouville problem. Similarly, the time-dependent functions $g_{1}(t)$
and $g_{2}(t)$ are assumed here to be at least once piece-wise differentiable.
Thus, they are continuous and their derivatives are piece-wise continuous on
$t\in[0,T]$. Standing on these hypotheses, a new method combining Fourier and
Laplace integral transformations is proposed here, in order to obtain
consistent analytical approximations in short time limits from the exact
operational solution to the boundary value problem (1)-(4). This implies that
all the involved functions are assumed here to be absolutely integrable with
respect to time and space variables, while the source term $f(x,t)$ is assumed
to be of exponential order relatively to the time variable $t$, as well as the
functions $g_{1}(t)$ and $g_{2}(t)$ and their respective derivatives.

##  3 Method of Fourier integral transform

In this section, in order to obtain the solution to the problem (1)-(4) under
a form of an integral, the Fourier integral transform (abbreviated FIT from
now on) and its inverse are used in relation to the space variable $x$. We
first recall that for any absolutely integrable function $\phi(x)$, i.e.,
$\int^{+\infty}_{-\infty}|\phi(x)|dx<\infty$, the FIT $\Phi$ is defined as:

$\Phi(\lambda)=\frac{1}{\sqrt{2\pi}}\int^{+\infty}_{-\infty}\phi(x)\exp(-i\lambda
x)dx,$ (5)

where $i^{2}=-1$, $\lambda\in\mathbb{R}$ and $x\in\mathbb{R}$. The inverse
Fourier transform of $\Phi$ is:

$\phi(x)=\frac{1}{\sqrt{2\pi}}\int^{+\infty}_{-\infty}\Phi(\lambda)\exp(ix\lambda)d\lambda.$
(6)

As a basic property, the FIT tends to 0 when $|\lambda|$ goes to $\infty$.

The homogeneous form of equation (1) reads:

$\frac{\partial u}{\partial t}=a^{2}\frac{\partial^{2}u}{\partial x^{2}}-bu.$
(7)

Let $u_{1}(x,t)$ be a solution to the boundary value problem formed by the
homogeneous equation (7) together with the initial condition (2), and the
boundary conditions (3) and (4). Assuming that the functions $u_{1}(x,t)$,
$\frac{\partial u_{1}}{\partial t}(x,t)$ and
$\frac{\partial^{2}u_{1}}{\partial x^{2}}(x,t)$ are absolutely integrable with
respect to the variables $x$ and $t$, they can be identified to their
extension by $0$ outside the rectangle $[l_{1},l_{2}]\times[0,T]$, without
lost of generality. The Fourier integral transform (FIT) with respect to the
space variable $x$ can be applied to $u_{1}(x,t)$ and will give:

$F(\lambda,t)=\frac{1}{\sqrt{2\pi}}\int^{+\infty}_{-\infty}u_{1}(x,t)\exp(-i\lambda
x)dx=\frac{1}{\sqrt{2\pi}}\int^{l_{2}}_{l_{1}}u_{1}(x,t)\exp(-i\lambda x)dx.$

The FIT applied to the transient term $\frac{\partial u_{1}}{\partial t}(x,t)$
leads to:

$A(\lambda,t)=\frac{1}{\sqrt{2\pi}}\int^{l_{2}}_{l_{1}}\frac{\partial
u_{1}}{\partial t}(x,t)\exp(-i\lambda
x)dx=\frac{1}{\sqrt{2\pi}}\frac{\partial}{\partial
t}\int^{l_{2}}_{l_{1}}u_{1}(x,t)\exp(-i\lambda x)dx=\frac{\partial F}{\partial
t}(\lambda,t)$

Similarly, the FIT of the diffusion term
$a^{2}\frac{\partial^{2}u_{1}}{\partial x^{2}}(x,t)$ can be written as:

$B(\lambda,t)=\frac{a^{2}}{\sqrt{2\pi}}\int^{l_{2}}_{l_{1}}\frac{\partial^{2}u_{1}}{\partial
x^{2}}(x,t)\exp(-i\lambda x)dx.$

By using two successive integration by parts over the finite interval
$[l_{1},l_{2}]$, $B(\lambda,t)$ can be expressed as:

$\begin{array}[]{ll}\displaystyle
B(\lambda,t)=-a^{2}\lambda^{2}F(\lambda,t)+\frac{a^{2}}{\sqrt{2\pi}}\left[\frac{\partial
u_{1}}{\partial x}(l_{2},t)\exp(-i\lambda l_{2})-\frac{\partial
u_{1}}{\partial x}(l_{1},t)\exp(-i\lambda l_{1})\right]\\\
\displaystyle+\frac{a^{2}}{\sqrt{2\pi}}\left[i\lambda
u_{1}(l_{2},t)\exp(-i\lambda l_{2})-i\lambda u_{1}(l_{1},t)\exp(-i\lambda
l_{1})\right].\end{array}$ (8)

The FIT of the linear term $-bu_{1}(x,t)$ is simply

$C(\lambda,t)=-bF(\lambda,t).$

Now, according to the application of the FIT to the terms of the homogeneous
equation, if $u_{1}(x,t)$, considered to be null outside
$[l_{1},l_{2}]\times[0,T]$ is solution of equation (7), then $F(\lambda,t)$ is
solution of the equation: $A(\lambda,t)-B(\lambda,t)-C(\lambda,t)=0$, i.e.,

$\begin{array}[]{ll}\displaystyle\frac{\partial F}{\partial
t}(\lambda,t)+(b+a^{2}\lambda^{2})F(\lambda,t)=\frac{a^{2}}{\sqrt{2\pi}}\left[\frac{\partial
u_{1}}{\partial x}(l_{2},t)\exp(-i\lambda l_{2})-\frac{\partial
u_{1}}{\partial x}(l_{1},t)\exp(-i\lambda l_{1})\right]\\\
\displaystyle+\frac{a^{2}}{\sqrt{2\pi}}\left[i\lambda
u_{1}(l_{2},t)\exp(-i\lambda l_{2})-i\lambda u_{1}(l_{1},t)\exp(-i\lambda
l_{1})\right];\end{array}$ (9)

for $\lambda\in\mathbb{R}$ and $0\leq t\leq T$.

In order to integrate equation (9), we multiply each member by
$\exp[(b+a^{2}\lambda^{2})t]$ and obtain:

$\begin{array}[]{ll}\displaystyle\frac{\partial}{\partial
t}(F(\lambda,t)\exp[(b+a^{2}\lambda^{2})t])=\\\
\displaystyle\frac{a^{2}}{\sqrt{2\pi}}\exp[(b+a^{2}\lambda^{2})t]\left[\frac{\partial
u_{1}}{\partial x}(l_{2},t)\exp(-i\lambda l_{2})-\frac{\partial
u_{1}}{\partial x}(l_{1},t)\exp(-i\lambda l_{1})\right]\\\
\displaystyle+\frac{a^{2}}{\sqrt{2\pi}}\exp[(b+a^{2}\lambda^{2})t]\left[i\lambda
u_{1}(l_{2},t)\exp(-i\lambda l_{2})-i\lambda u_{1}(l_{1},t)\exp(-i\lambda
l_{1})\right].\end{array}$ (10)

Proceeding now by integration of equation (10) relatively to the time variable
from $\eta=0$ to $\eta=t\leq T$, we obtain:

$\begin{array}[]{ll}\displaystyle
F(\lambda,t)\exp[(b+a^{2}\lambda^{2})t]-F(\lambda,0)=\\\
\displaystyle\frac{a^{2}}{\sqrt{2\pi}}\int_{0}^{t}\left[\frac{\partial
u_{1}}{\partial x}(l_{2},\eta)\exp(-i\lambda l_{2})-\frac{\partial
u_{1}}{\partial x}(l_{1},\eta)\exp(-i\lambda
l_{1})\right]\exp[(b+a^{2}\lambda^{2})\eta]d\eta\\\
\displaystyle+\frac{a^{2}}{\sqrt{2\pi}}\int_{0}^{t}\left[i\lambda
u_{1}(l_{2},\eta)\exp(-i\lambda l_{2})-i\lambda
u_{1}(l_{1},\eta)\exp(-i\lambda
l_{1})\right]\exp[(b+a^{2}\lambda^{2})\eta]d\eta.\end{array}$ (11)

Due to the initial condition (2), $F(\lambda,0)=F(\lambda,t=0)$ can be
calculated as

$F(\lambda,0)=\frac{1}{\sqrt{2\pi}}\int^{l_{2}}_{l_{1}}\varphi(\xi)\exp(-i\lambda\xi)d\xi,$
(12)

where the dummy variable is replaced by the integrating variable $\xi$ in
order to avoid confusion. Equation (11) can then be rewritten as:

$\begin{array}[]{ll}\displaystyle F(\lambda,t)=\\\
\displaystyle\frac{a^{2}}{\sqrt{2\pi}}\int_{0}^{t}\left[\frac{\partial
u_{1}}{\partial x}(l_{2},\eta)\exp(-i\lambda l_{2})-\frac{\partial
u_{1}}{\partial x}(l_{1},\eta)\exp(-i\lambda
l_{1})\right]\exp[-(b+a^{2}\lambda^{2})(t-\eta)]d\eta\\\
\displaystyle+\frac{a^{2}}{\sqrt{2\pi}}\int_{0}^{t}\left[i\lambda
u_{1}(l_{2},\eta)\exp(-i\lambda l_{2})-i\lambda
u_{1}(l_{1},\eta)\exp(-i\lambda
l_{1})\right]\exp[-(b+a^{2}\lambda^{2})(t-\eta)]d\eta\\\
\displaystyle+\frac{1}{\sqrt{2\pi}}\exp[-(b+a^{2}\lambda^{2})t]\int^{l_{2}}_{l_{1}}\varphi(\xi)\exp(-i\lambda\xi)d\xi.\end{array}$
(13)

In order to obtain $u_{1}(x,t)$, the inversion formula (6) needs to be applied
to the function $F(\lambda,t)$ expressed by equation (13), the Fourier
variable $\lambda$ running from $-\infty$ to $+\infty$. First, due to
properties of the convolution of two functions, the three terms of the second
member can be respectively denoted as:

$F_{1}(\lambda,t)=\frac{a^{2}}{\sqrt{2\pi}}\int_{0}^{t}\left[\frac{\partial
u_{1}}{\partial x}(l_{2},t-\eta)\exp(-i\lambda l_{2})-\frac{\partial
u_{1}}{\partial x}(l_{1},t-\eta)\exp(-i\lambda
l_{1})\right]\exp[-(b+a^{2}\lambda^{2})\eta]d\eta,$
$F_{2}(\lambda,t)=\frac{a^{2}}{\sqrt{2\pi}}\int_{0}^{t}\left[i\lambda
u_{1}(l_{2},t-\eta)\exp(-i\lambda l_{2})-i\lambda
u_{1}(l_{1},t-\eta)\exp(-i\lambda
l_{1})\right]\exp[-(b+a^{2}\lambda^{2})\eta]d\eta$

and

$F_{3}(\lambda,t)=\frac{1}{\sqrt{2\pi}}\exp[-(b+a^{2}\lambda^{2})t]\int^{l_{2}}_{l_{1}}\varphi(\xi)\exp(-i\lambda\xi)d\xi.$

Changing the order of integration due to the convergence of the integrals
involved, the inverse $I_{1}$ of the first term $F_{1}$ is calculated as:

$\begin{array}[]{ll}I_{1}(x,t)=\displaystyle\frac{a^{2}}{{2\pi}}\int_{0}^{t}\frac{\partial
u_{1}}{\partial x}(l_{2},t-\eta)\int_{-\infty}^{\infty}\exp(-i\lambda
l_{2})\exp(i\lambda x)\exp[-(b+a^{2}\lambda^{2})\eta]d\lambda d\eta\\\
-\displaystyle\frac{a^{2}}{{2\pi}}\int_{0}^{t}\frac{\partial u_{1}}{\partial
x}(l_{1},t-\eta)\int_{-\infty}^{\infty}\exp(-i\lambda l_{1})\exp(i\lambda
x)\exp[-(b+a^{2}\lambda^{2})\eta]d\lambda d\eta.\end{array}$

By means of computations $I_{1}$ is reduced to:

$\begin{array}[]{ll}I_{1}(x,t)=\displaystyle\frac{a}{2\sqrt{\pi}}\int_{0}^{t}\frac{\partial
u_{1}}{\partial
x}(l_{2},t-\eta)\left[\exp\left(-\frac{(l_{2}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{\sqrt{\eta}}d\eta\\\
-\displaystyle\frac{a}{2\sqrt{\pi}}\int_{0}^{t}\frac{\partial u_{1}}{\partial
x}(l_{1},t-\eta)\left[\exp\left(-\frac{(l_{1}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{\sqrt{\eta}}d\eta.\end{array}$
(14)

Similarly, $I_{2}$ is obtained as:

$\begin{array}[]{ll}I_{2}(x,t)=\displaystyle\frac{a^{2}}{{2\pi}}\int_{0}^{t}u_{1}(l_{2},t-\eta)\int_{-\infty}^{\infty}i\lambda\exp(-i\lambda
l_{2})\exp(i\lambda x)\exp[-(b+a^{2}\lambda^{2})\eta]d\lambda d\eta\\\
-\displaystyle\frac{a^{2}}{{2\pi}}\int_{0}^{t}u_{1}(l_{1},t-\eta)\int_{-\infty}^{\infty}i\lambda\exp(-i\lambda
l_{1})\exp(i\lambda x)\exp[-(b+a^{2}\lambda^{2})\eta]d\lambda
d\eta,\end{array}$

i.e.,

$\begin{array}[]{ll}I_{2}(x,t)=\displaystyle\frac{1}{4a\sqrt{\pi}}\int_{0}^{t}u_{1}(l_{2},t-\eta)\left[(l_{2}-x)\exp\left(-\frac{(l_{2}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{{\eta}^{3/2}}d\eta\\\
-\displaystyle\frac{1}{4a\sqrt{\pi}}\int_{0}^{t}u_{1}(l_{1},t-\eta)\left[(l_{1}-x)\exp\left(-\frac{(l_{1}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{{\eta}^{3/2}}d\eta.\end{array}$
(15)

And, the inverse $I_{3}$ of $F_{3}$ is calculated as:

$I_{3}(x,t)=\displaystyle\frac{\exp(-bt)}{2a\sqrt{\pi
t}}\int_{l_{1}}^{l_{2}}\varphi(\xi)\left[\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}t}\right)\right]d\xi.$
(16)

From the calculations above, the expression of $u_{1}$ is deduced as
$u_{1}(x,t)=I_{1}(x,t)+I_{2}(x,t)+I_{3}(x,t)$, that is:

$\begin{array}[]{ll}u_{1}(x,t)=\displaystyle\frac{a}{2\sqrt{\pi}}\int_{0}^{t}\frac{\partial
u_{1}}{\partial
x}(l_{2},t-\eta)\left[\exp\left(-\frac{(l_{2}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{\sqrt{\eta}}d\eta\\\
-\displaystyle\frac{a}{2\sqrt{\pi}}\int_{0}^{t}\frac{\partial u_{1}}{\partial
x}(l_{1},t-\eta)\left[\exp\left(-\frac{(l_{1}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{\sqrt{\eta}}d\eta\\\
+\displaystyle\frac{1}{4a\sqrt{\pi}}\int_{0}^{t}u_{1}(l_{2},t-\eta)\left[(l_{2}-x)\exp\left(-\frac{(l_{2}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{{\eta}^{3/2}}d\eta\\\
-\displaystyle\frac{1}{4a\sqrt{\pi}}\int_{0}^{t}u_{1}(l_{1},t-\eta)\left[(l_{1}-x)\exp\left(-\frac{(l_{1}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{{\eta}^{3/2}}d\eta\\\
+\displaystyle\frac{\exp(-bt)}{2a\sqrt{\pi
t}}\int_{l_{1}}^{l_{2}}\varphi(\xi)\left[\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}t}\right)\right]d\xi.\end{array}$
(17)

It can be verified that $u_{1}$, as expressed by the integral form (17), is a
solution to the homogeneous equation (7). Moreover, the initial condition (2)
is satisfied by the solution (17) since, when $t\rightarrow 0$, $I_{1}(x,t)$
and $I_{2}(x,t)$ vanish and $\displaystyle\lim_{t\rightarrow 0}u_{1}(x,t)$
reduces to $\displaystyle\lim_{t\rightarrow 0}I_{3}(x,t)$. Now, if we put

$G(x,\xi,t)=\frac{\exp(-bt)}{2a\sqrt{\pi
t}}\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}t}\right),$

then, due to the property of the normal or Gauss probability density:

$\displaystyle\lim_{t\rightarrow
0}\int_{-\infty}^{\infty}G(x,\xi,t)d\xi=\displaystyle\lim_{t\rightarrow
0}\exp(-bt)=1,$

and $G(x,\xi,t)\rightarrow 0$ as $t\rightarrow 0$ at all points
$(x,\xi)\in\mathbb{R}^{2}$, with the exception of the diagonal $x=\xi$ where
it becomes infinitely large. Therefore, if $G(x,\xi,t)$ is an analogue of a
Green’s function and

$\displaystyle\lim_{t\rightarrow 0}G(x,\xi,t)=\delta(x-\xi),$

where $\delta(x-\xi)$ is the Dirac delta function. If
$\mathbb{I}_{[l_{1},l_{2}]}$ denotes the indicator function of interval
$[l_{1},l_{2}]$, one has:

$\begin{array}[]{ll}\displaystyle\lim_{t\rightarrow
0}u_{1}(x,t)=\displaystyle\lim_{t\rightarrow
0}I_{3}(x,t)=\displaystyle\lim_{t\rightarrow
0}\int_{l_{1}}^{l_{2}}\varphi(\xi)G(x,\xi,t)d\xi\\\
=\displaystyle\lim_{t\rightarrow
0}\int_{-\infty}^{\infty}\mathbb{I}_{[l_{1},l_{2}]}(\xi)\varphi(\xi)G(x,\xi,t)d\xi=\int_{-\infty}^{\infty}\mathbb{I}_{[l_{1},l_{2}]}(\xi)\varphi(\xi)\delta(x-\xi)d\xi=\varphi(x),\end{array}$
(18)

and the initial condition (2) is satisfied by $u_{1}$. In fact, it has just
been proven that the function $I_{3}$ satisfied the non-zero initial condition
(2). Moreover, it can also be verified that $I_{3}$ is a solution to the
homogeneous equation (7). According to the Duhamel’s principle, a solution of
the nonhomogeneous equation (1) with zero initial condition can be written as:

$u_{2}(x,t)=\int_{0}^{t}d\theta\int_{l_{1}}^{l_{2}}\,G(x,\xi,t-\theta)f(\xi,\theta)d\xi.$

Finally, by the superposition principle, a solution $u$ of the nonhomogeneous
equation (1) together with the non-zero initial condition (2) can be expressed
as $u=u_{1}+u_{2}$, namely:

$\begin{array}[]{ll}u(x,t)=\displaystyle\frac{a}{2\sqrt{\pi}}\int_{0}^{t}\frac{\partial
u}{\partial
x}(l_{2},t-\eta)\left[\exp\left(-\frac{(l_{2}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{\sqrt{\eta}}d\eta\\\
-\displaystyle\frac{a}{2\sqrt{\pi}}\int_{0}^{t}\frac{\partial u}{\partial
x}(l_{1},t-\eta)\left[\exp\left(-\frac{(l_{1}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{\sqrt{\eta}}d\eta\\\
+\displaystyle\frac{1}{4a\sqrt{\pi}}\int_{0}^{t}u(l_{2},t-\eta)\left[(l_{2}-x)\exp\left(-\frac{(l_{2}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{{\eta}^{3/2}}d\eta\\\
-\displaystyle\frac{1}{4a\sqrt{\pi}}\int_{0}^{t}u(l_{1},t-\eta)\left[(l_{1}-x)\exp\left(-\frac{(l_{1}-x)^{2}}{4a^{2}\eta}\right)\right]\frac{\exp(-b\eta)}{{\eta}^{3/2}}d\eta\\\
+r(x,t)\end{array}$ (19)

where:

$\begin{array}[]{ll}r(x,t)=\displaystyle\frac{\exp(-bt)}{2a\sqrt{\pi
t}}\int_{l_{1}}^{l_{2}}\varphi(\xi)\left[\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}t}\right)\right]d\xi\\\
+\displaystyle\frac{1}{2a\sqrt{\pi}}\int_{0}^{t}d\theta\int_{l_{1}}^{l_{2}}\frac{\exp(-b(t-\theta))}{\sqrt{(t-\theta)}}\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}(t-\theta)}\right)f(\xi,\theta)d\xi\end{array}$
(20)

The FIT method has allowed to determine an expression in integral form of the
solution of the equation (1) under the initial condition (2). But, the
solution (19) depends of the values of the function $u$ and of its derivative
at the boundaries. The boundary conditions (3) and (4) will be taken into
account via the Laplace integral transform (abbreviated LIT).

##  4 Exact solution in the Laplace domain

If $f(t)$ is a function defined in $t\geq 0$, then its unilateral Laplace
integral transform (LIT) is given in the complex $p$-plane by (see [14]):

$F(p)={\mathcal{L}}\\{f(t)\\}=\int_{0}^{\infty}f(t)\mathrm{e}^{-pt}dt,$ (21)

provided that $f(t)$ be of exponential order, that is, there are constants $C$
and $\sigma$ so that $|f(t)|<C\mathrm{e}^{\sigma t}$, when t is sufficiently
large. The inversion, from the Laplace domain $p$ to the time domain $t$ is
given by the complex integral,

$f(t)={\mathcal{L}}^{-1}\\{F(p)\\}=\frac{1}{2\pi
i}\int_{\gamma-i\infty}^{\gamma+i\infty}F(p)\mathrm{e}^{pt}dp,$ (22)

where $\gamma>\sigma$ is chosen so that $F(p)$ converges absolutely on the
real part of $p$ line $\Re(p)=\gamma$, and $F(p)$ is analytic at the right of
this line. Also, the LIT and its inverse can be obtained for many usual
functions by using tables of transforms as in [24]. In the case when $p$ is
real as considered in the sequel of this paper, the inequality $p\geq\sigma$
needs to be satisfied. Relatively to the LIT, an important property is the
convolution theorem (see [8]):
Let $f(t)$ and $g(t)$ be functions defined in $t\geq 0$. If
${\mathcal{L}}\\{f(t)\\}=F(p)$ and ${\mathcal{L}}\\{g(t)\\}=G(p)$, then

${\mathcal{L}}\\{f(t)\ast
g(t)\\}={\mathcal{L}}\\{f(t)\\}{\mathcal{L}}\\{g(t)\\}=F(p)G(p)$

where $f(t)\ast g(t)$ is called the convolution of $f(t)$ and $g(t)$ and is
defined by the integral

$f(t)\ast g(t)=\int^{t}_{0}f(t-\eta)g(\eta)d\eta.$

In order to apply the LIT, all the time dependent functions involved in the
problem (1)-(4) are assumed to be originals, that is, the transforms of these
functions exist. Thus, the LIT of temperature distribution $u(x,t)$, source
term $f(x,t)$, boundary functions $g_{1}(t)$ and $g_{2}(t)$ are respectively
denoted in the Laplace domain by $U(x,p)$, $F(x,p)$, $G_{1}(p)$ and
$G_{2}(p)$. Under these hypotheses, the analogue of the initial boundary value
problem (1)-(4) in the Laplace domain can be written in the form of an
ordinary differential equation as follows:

$-a^{2}\frac{d^{2}U}{dx^{2}}(x,p)+(b+p)U(x,p)=F(x,p)+\varphi(x),$ (23)

since the Laplace transform of the time derivative is equal to the product of
the transform by the operator p minus the value of the function at the initial
time instant:

${\mathcal{L}}\left\\{\frac{\partial u}{\partial
t}(x,t)\right\\}=pU(x,p)-\varphi(x).$ (24)

Equation (23) is subjected to the boundary conditions:

$\alpha_{1}U(l_{1},p)+\beta_{1}\frac{dU}{dx}(l_{1},p)=G_{1}(p),\quad\alpha_{1},\beta_{1}\in\mathbb{R},\quad\alpha_{1}^{2}+\beta_{1}^{2}\neq
0,$ (25)
$\alpha_{2}U(l_{2},p)+\beta_{2}\frac{dU}{dx}(l_{2},p)=G_{2}(p),\quad\alpha_{2},\beta_{2}\in\mathbb{R},\quad\alpha_{2}^{2}+\beta_{2}^{2}\neq
0.$ (26)

Due to the convolution theorem, the solution expressed by equation (19) can
also be transformed by the LIT to be written as:

$\begin{array}[]{ll}U(x,p)=\\\
\displaystyle\frac{a}{2\sqrt{b+p}}\left[U_{x}(l_{2},p)\exp\left(\frac{-(l_{2}-x)\sqrt{b+p}}{a}\right)-U_{x}(l_{1},p)\exp\left(\frac{-(x-l_{1})\sqrt{b+p}}{a}\right)\right]\\\
+\displaystyle\frac{1}{2}\left[U(l_{2},p)\exp\left(\frac{-(l_{2}-x)\sqrt{b+p}}{a}\right)+U(l_{1},p)\exp\left(\frac{-(x-l_{1})\sqrt{b+p}}{a}\right)\right]\\\
+\displaystyle R(x,p),\end{array}$ (27)

where $U_{x}$ denotes the derivative $\frac{dU}{dx}$ and $R(x,p)$ stands for
the Laplace transform of the remaining term $r(x,t)$ expressed by the equation
(20), namely, $R(x,p)={\mathcal{L}}\\{r(x,t)\\}$. Substituting respectively
$x=l_{1}$ and then $x=l_{2}$ in the above expression of $U(x,p)$, we have the
following equations:

$\begin{array}[]{ll}\displaystyle\frac{1}{2}U(l_{1},p)+\frac{a}{2\sqrt{b+p}}U_{x}(l_{1},p)-\frac{1}{2}\exp\left(\frac{-(l_{2}-l_{1})\sqrt{b+p}}{a}\right)U(l_{2},p)\\\
\displaystyle-\frac{a}{2\sqrt{b+p}}\exp\left(\frac{-(l_{2}-l_{1})\sqrt{b+p}}{a}\right)U_{x}(l_{2},p)=R(l_{1},p),\end{array}$
(28)

and

$\begin{array}[]{ll}\displaystyle-\frac{1}{2}\exp\left(\frac{-(l_{2}-l_{1})\sqrt{b+p}}{a}\right)U(l_{1},p)+\frac{a}{2\sqrt{b+p}}\exp\left(\frac{-(l_{2}-l_{1})\sqrt{b+p}}{a}\right)U_{x}(l_{1},p)\\\
\displaystyle+\frac{1}{2}U(l_{2},p)-\frac{a}{2\sqrt{b+p}}U_{x}(l_{2},p)=R(l_{2},p),\end{array}$
(29)

Now, equations (25), (26), (28) and (29) form a system ($\mathcal{S}$) of four
linear equations with four unknown functions that are $U(l_{1},p)$,
$U_{x}(l_{1},p)$, $U(l_{2},p)$ and $U(l_{2},p)$. Therefore, the function
$U(x,p)$ is solution of the analogue problem (23)-(26) in the Laplace domain,
or equivalently $u(x,t)$ is solution of the problem (1)-(4) in the time
domain, if and only if, the system ($\mathcal{S}$) admits a unique solution.
The determinant of the system is calculated as:

$\mathrm{det}(\mathcal{S})=\left|\begin{array}[]{cccc}\displaystyle\alpha_{1}&\displaystyle\beta_{1}&\displaystyle
0&\displaystyle 0\\\ \displaystyle
0&0&\displaystyle\alpha_{2}&\displaystyle\beta_{2}\\\
\displaystyle\frac{1}{2}&\displaystyle\frac{a}{2\sqrt{b+p}}&\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle\frac{-a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}\\\
\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle\frac{a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}&\displaystyle\frac{1}{2}&\displaystyle\frac{-a}{2\sqrt{b+p}}\end{array}\right|;$
(30)

where the function $\chi(x,p)$ is defined as
$\chi(x,p)=\exp\left(\frac{-x\sqrt{b+p}}{a}\right)$. This determinant is
reduced by computations into:

$\begin{array}[]{cccc}\displaystyle\det(\mathcal{S})=-\frac{1}{4}(a^{2}\alpha_{1}\alpha_{2}-a\alpha_{1}\beta_{2}+a\alpha_{2}\beta_{1})\frac{\chi(2(l_{2}-l_{1}),p)}{b+p}+\frac{1}{4}\beta_{1}\beta_{2}\chi(2(l_{2}-l_{1}),p)\\\
\displaystyle+\frac{1}{4}\frac{a^{2}\alpha_{1}\alpha_{2}+a\alpha_{1}\beta_{2}-a\alpha_{2}\beta_{1}}{b+p}-\frac{1}{4}\beta_{1}\beta_{2}.\end{array}$
(31)

Thus, the determinant of the system is null if and only if all the
coefficients of the functions appearing in equation (31) are null. This will
imply:

$\left\\{\begin{array}[]{cccc}a^{2}\alpha_{1}\alpha_{2}-a\alpha_{1}\beta_{2}+a\alpha_{2}\beta_{1}=0,\\\
a^{2}\alpha_{1}\alpha_{2}+a\alpha_{1}\beta_{2}-a\alpha_{2}\beta_{1}=0,\\\
\beta_{1}\beta_{2}=0.\end{array}\right.$ (32)

The system (32) can be split into two systems since it is equivalent to:

$\left\\{\begin{array}[]{cccc}a^{2}\alpha_{1}\alpha_{2}-a\alpha_{1}\beta_{2}=0,\\\
a^{2}\alpha_{1}\alpha_{2}+a\alpha_{1}\beta_{2}=0,\\\
\beta_{1}=0,\end{array}\right.$ (33)

or

$\left\\{\begin{array}[]{cccc}a^{2}\alpha_{1}\alpha_{2}+a\alpha_{2}\beta_{1}=0,\\\
a^{2}\alpha_{1}\alpha_{2}-a\alpha_{2}\beta_{1}=0,\\\
\beta_{2}=0.\end{array}\right.$ (34)

Since $a\neq 0$, by adding and then subtracting the two first equations, the
sub-system (33) can be shown to be equivalent to:

$\left\\{\begin{array}[]{cccc}\alpha_{1}\alpha_{2}=0,\\\
\alpha_{1}\beta_{2}=0,\\\ \beta_{1}=0.\end{array}\right.$ (35)

The system (35) leads to a contradiction with the hypotheses on the
coefficients $\alpha_{1},\>\beta_{1},\>\alpha_{2}$ and $\beta_{2}$, since it
implies that $\alpha_{1}=\beta_{1}=0$ or $\alpha_{2}=\beta_{2}=0$. The same
contradiction is reached when trying to solve the system (34). So,
$\det(\mathcal{S})\neq 0$ in all cases. Consequently, the analogue boundary
value problem (23)-(26) in the Laplace domain, or equivalently the problem
(1)-(4) in the time domain, admits a unique solution whenever
$\alpha_{1}^{2}+\beta_{1}^{2}\neq 0$ and $\alpha_{2}^{2}+\beta_{2}^{2}\neq 0$.
Now, the exact solution in the Laplace domain of the problem (1)-(4) is
$U(x,p)$ given by equation (27), with the functions $U(l_{1},p)$,
$U(l_{2},p)$, $U_{x}(l_{1},p)$ and $U_{x}(l_{2},p)$ expressed by using
determinants as:

$U(l_{1},p)=\frac{\left|\begin{array}[]{cccc}\displaystyle
G_{1}(p)&\displaystyle\beta_{1}&\displaystyle 0&\displaystyle 0\\\
\displaystyle G_{2}(p)&0&\displaystyle\alpha_{2}&\displaystyle\beta_{2}\\\
\displaystyle
R(l_{1},p)&\displaystyle\frac{a}{2\sqrt{b+p}}&\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle\frac{-a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}\\\
\displaystyle
R(l_{2},p)&\displaystyle\frac{a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}&\displaystyle\frac{1}{2}&\displaystyle\frac{-a}{2\sqrt{b+p}}\end{array}\right|}{\mathrm{det}(\mathcal{S})};$
(36)
$U_{x}(l_{1},p)=\frac{\left|\begin{array}[]{cccc}\displaystyle\alpha_{1}&\displaystyle
G_{1}(p)&\displaystyle 0&\displaystyle 0\\\ \displaystyle
0&G_{2}(p)&\displaystyle\alpha_{2}&\displaystyle\beta_{2}\\\
\displaystyle\frac{1}{2}&\displaystyle
R(l_{1},p)&\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle\frac{-a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}\\\
\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle
R(l_{2},p)&\displaystyle\frac{1}{2}&\displaystyle\frac{-a}{2\sqrt{b+p}}\end{array}\right|}{\mathrm{det}(\mathcal{S})};$
(37)
$U(l_{2},p)=\frac{\left|\begin{array}[]{cccc}\displaystyle\alpha_{1}&\displaystyle\beta_{1}&\displaystyle
G_{1}(p)&\displaystyle 0\\\ \displaystyle 0&0&\displaystyle
G_{2}(p)&\displaystyle\beta_{2}\\\
\displaystyle\frac{1}{2}&\displaystyle\frac{a}{2\sqrt{b+p}}&\displaystyle
R(l_{1},p)&\displaystyle\frac{-a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}\\\
\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle\frac{a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}&\displaystyle
R(l_{2},p)&\displaystyle\frac{-a}{2\sqrt{b+p}}\end{array}\right|}{\mathrm{det}(\mathcal{S})};$
(38)

and

$U_{x}(l_{x},p)=\frac{\left|\begin{array}[]{cccc}\displaystyle\alpha_{1}&\displaystyle\beta_{1}&\displaystyle
0&\displaystyle G_{1}(p)\\\ \displaystyle
0&0&\displaystyle\alpha_{2}&\displaystyle G_{2}(p)\\\
\displaystyle\frac{1}{2}&\displaystyle\frac{a}{2\sqrt{b+p}}&\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle
R(l_{1},p)\\\
\displaystyle\frac{-\chi(l_{2}-l_{1},p)}{2}&\displaystyle\frac{a\chi(l_{2}-l_{1},p)}{2\sqrt{b+p}}&\displaystyle\frac{1}{2}&\displaystyle
R(l_{2},p)\end{array}\right|}{\mathrm{det}(\mathcal{S})}.$ (39)

In brief, the exact solution in the Laplace domain of the boundary value
problem (1)-(4) is expressed in a unified way by the function $U(x,p)$ given
by equation (27) together with those given in equations (36)-(39). Exact
series solutions in the time domain of such linear boundary value problems are
well established by using the Fourier decomposition method as in [2] for
example. Thus, the exact Laplace transform of those series solutions performed
via the Sturm-Liouville theory can be recovered by the expression of the
function $U(x,p)$.

##  5 Analytical approximations in short time limits

In this section, our interest is to find approximate analytical solutions to
the problem (1)-(4) at the earliest times of the process. We assume a
subdivision of the time interval $[0,T]$, such that the magnitude of the
dimensionless time step $\Delta t$ is sufficiently small, for example $\Delta
t\leq 10^{-2}$, as in many schemes in numerical analysis (see [25] and [15]
among others). The goal of this section is to derive from the exact
operational solution (27), approximate analytical solutions valid during the
first time step of the reaction–diffusion process, namely for $t\in[0,\Delta
t]$.

The limiting case of a short time ($\Delta t$ tending to 0) corresponds to a
very large value of the Laplace domain variable ($p$ tending to $+\infty$). In
order to minimize calculations, some prior simplifications can be done on the
determinant of the system $\mathrm{det}(\mathcal{S})$, before deriving the
asymptotic expansions of the solution and the related truncated expansions in
the time domain. For $\eta>0$, the inverse Laplace transform of $\chi(\eta,p)$
is written as:

${\mathcal{L}}^{-1}\\{\chi(\eta,p)\\}={\mathcal{L}}^{-1}\left\\{\exp\left(\frac{-\eta\sqrt{b+p}}{a}\right)\right\\}=\frac{1}{2}\frac{\eta\exp\left(-bt-\frac{\eta^{2}}{4a^{2}t}\right)}{a\sqrt{\pi}t^{3/2}}=\Gamma(t).$
(40)

When $p\rightarrow\infty$, $\chi(\eta,p)$ is negligible compared to
${1}/{p^{n}}$, for all positive integers $n\geq 1$, whereas when $t\rightarrow
0$, the inverse Laplace transform $\Gamma(t)$ is negligible compared to
$t^{n}$. Indeed,
$\displaystyle\lim_{p\rightarrow\infty}p^{n}\chi(\eta,p)=\displaystyle\lim_{t\rightarrow
0}\Gamma(t)/t^{n}=0$, implying $\chi(\eta,p)=o({1}/{p^{n}})$ at $p=\infty$ and
$\Gamma(t)=o(t^{n})$ at $t=0$, where the Little-o is the asymptotic notation.
Again, these relations are valid for arbitrary order $n\geq 1$. When
$p\rightarrow\infty$, the term $\chi(l_{2}-l_{1},p)$ and its asymptotic
expansions are negligible in the Laplace domain as well as their inverses in
the time domain when $t\rightarrow 0$. Therefore, when $p$ is sufficiently
large, the asymptotic solution can be obtained using the following determinant
for the system formed by equations (25), (26), (28) and (29):

$\mathrm{det}(\mathcal{S}^{a})=\left|\begin{array}[]{cccc}\displaystyle\alpha_{1}&\displaystyle\beta_{1}&\displaystyle
0&\displaystyle 0\\\ \displaystyle
0&0&\displaystyle\alpha_{2}&\displaystyle\beta_{2}\\\
\displaystyle\frac{1}{2}&\displaystyle\frac{a}{2\sqrt{b+p}}&\displaystyle
0&\displaystyle 0\\\ \displaystyle 0&\displaystyle
0&\displaystyle\frac{1}{2}&\displaystyle\frac{-a}{2\sqrt{b+p}}\end{array}\right|.$
(41)

Based on the same assumption that the function $\chi(\eta,p)$ is negligible
for all $\eta>0$ whenever $p$ is sufficiently large, an asymptotic expansion
of the exact solution $U(x,p)$, valid at any order $n\geq 1$, can be deduced
for $l_{1}<x<l_{2}$, from equation (27) as follows:

$U^{a}(x,p)=R(x,p)+o\left(\frac{1}{p^{n}}\right).$ (42)

Equation (42) corresponds for $0<t<\Delta t$ to the following truncated
expansion of the solution $u(x,t)$ in the time domain:

$u^{a}(x,t)=r(x,t)+o\left((\Delta t)^{n}\right);$ (43)

where $r(x,t)$ is given by equation (20) and $o\left((\Delta t)^{n}\right)$ is
an upper bound of the truncation error $o\left(t^{n}\right)$. But, at
$x=l_{1}$ and $x=l_{2}$, some refined asymptotic expansions in the Laplace
domain $U^{a}(l_{1},p)$, $U^{a}(l_{2},p)$, $U^{a}_{x}(l_{1},p)$ and
$U^{a}_{x}(l_{2},p)$ can be obtained for the boundary-related functions
$U(l_{1},p)$, $U(l_{2},p)$, $U_{x}(l_{1},p)$ and $U_{x}(l_{2},p)$. We
distinguish four different cases with respect to the values of the two
coefficients $\beta_{1}$ and $\beta_{2}$. For each case, the corresponding
analytical approximations in the time domain, namely $u^{a}(l_{1},t)$,
$u^{a}(l_{2},t)$, $u^{a}_{x}(l_{1},t)$ and $u^{a}_{x}(l_{2},t)$, are also
given for $t\in[0,\Delta t]$. Although this can be improved, the asymptotic
expansions at second order are retained here, while the corresponding
analytical approximations in the time domain are given at the first order.

  * •

Case $\beta_{1}\beta_{2}\neq 0$

    1. 1.

Solving the system by using the reduced determinant
$\mathrm{det}(\mathcal{S}^{a})$ expressed above in formula (41), leads to:

$\displaystyle
U(l_{1},p)=-\frac{a}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}\,G_{1}(p)+\frac{2\sqrt{b+p}\beta_{1}}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}\,R(l_{1},p),$

and if the Big-O denotes the asymptotic notation, it can be written:

$\displaystyle-\frac{a}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}=U^{1l_{1}}(p)+O\left(\frac{1}{p^{2}}\right);$

and

$\displaystyle\frac{2\sqrt{b+p}\beta_{1}}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}=2+U^{2l_{1}}(p)+O\left(\frac{1}{p^{2}}\right);$

where

$U^{1l_{1}}(p)=\displaystyle-\frac{a}{\beta_{1}}\sqrt{\frac{1}{p}}-{\frac{{a}^{2}\alpha_{1}}{{\beta_{1}}^{2}p}}-{\frac{1}{\beta_{1}}\left(-\frac{1}{2}\,ab+{\frac{{a}^{3}{\alpha_{1}}^{2}}{{\beta_{1}}^{2}}}\right)\left(\frac{1}{p}\right)^{3/2}};$

and

$U^{2l_{1}}(p)=\displaystyle
2\,{\frac{a\alpha_{1}}{\beta_{1}}}\sqrt{\frac{1}{p}}+2\,\frac{{a}^{2}{\alpha_{1}}^{2}}{{\beta_{1}}^{2}p}+\,{\frac{2}{\beta_{1}}\left(-\frac{1}{2}\,a\alpha_{1}b+{\frac{{a}^{3}{\alpha_{1}}^{3}}{{\beta_{1}}^{2}}}\right)\left(\frac{1}{p}\right)^{3/2}}.$

Now, due to the properties of the LIT, $G_{1}(p)$ and $R(l_{1},p)$ are bounded
functions in the Laplace domain for $p>0$, and an asymptotic expansion of the
solution $U(l_{1},p)$ can be written as:

$\displaystyle
U^{a}(l_{1},p)=U^{1l_{1}}(p)G_{1}(p)+U^{2l_{1}}(p)R(l_{1},p)+2R(l_{1},p)+O\left(\frac{1}{p^{2}}\right).$

The corresponding analytical approximation in the time domain during a short
time step $t\in[0,\Delta t]$ reads

$u^{a}(l_{1},t)=u^{1l_{1}}(t)\ast g_{1}(t)+u^{2l_{1}}(t)\ast
r(l_{1},t)+2r(l_{1},t)+O\left(\Delta t\right),$

where $O\left(\Delta t\right)$ instead of $O\left(t\right)$ is an upper bound
of the truncation error of the approximation, $\ast$ denotes the convolution
product, $r(l_{1},t)$ is calculated through equation (20), $g_{1}(t)$ is the
known function related to the boundary $x=l_{1}$,

$u^{1l_{1}}(t)={\mathcal{L}}^{-1}\\{U^{1l_{1}}(p)\\}=-{\frac{{a}^{2}\alpha_{1}}{{\beta_{1}}^{2}}}-{\frac{a}{\sqrt{\pi\,t}\beta_{1}}}+{\frac{\left(-2\,{a}^{2}{\alpha_{1}}^{2}+b{\beta_{1}}^{2}\right)a}{{\beta_{1}}^{3}}\sqrt{{\frac{t}{\pi}}}};$

and

$u^{2l_{1}}(t)={\mathcal{L}}^{-1}\\{U^{2l_{1}}(p)\\}=2\,{\frac{a\alpha_{1}}{\sqrt{\pi\,t}\beta_{1}}}+2\,{\frac{{a}^{2}{\alpha_{1}}^{2}}{{\beta_{1}}^{2}}}+2\,{\frac{a\alpha_{1}\left(2\,{a}^{2}{\alpha_{1}}^{2}-b{\beta_{1}}^{2}\right)}{{\beta_{1}}^{3}}\sqrt{{\frac{t}{\pi}}}}.$

    2. 2.

For $U(l_{2},p)$ and $u(l_{2},t)$, similar calculations as above give:

$U(l_{2},p)={\frac{G_{2}\left(p\right)a}{\beta_{2}\sqrt{b+p}+\alpha_{{2}}a}}+2\,{\frac{R\left(l_{2},p\right)\sqrt{b+p}\beta_{2}}{\beta_{2}\sqrt{b+p}+\alpha_{2}a}};$

and

$\displaystyle
U^{a}(l_{2},p)=U^{1l_{2}}(p)G_{2}(p)+U^{2l_{2}}(p)R(l_{2},p)+2R(l_{2},p)+O\left(\frac{1}{p^{2}}\right);$

where

$U^{1l_{2}}(p)=\displaystyle{\frac{a}{\beta_{2}}}\sqrt{\frac{1}{p}}-{\frac{{a}^{2}\alpha_{2}}{{\beta_{2}}^{2}p}}+{\frac{1}{\beta_{{2}}}\left(-\frac{1}{2}\,ab+{\frac{{a}^{3}{\alpha_{2}}^{2}}{{\beta_{2}}^{2}}}\right)\left(\frac{1}{p}\right)^{3/2}};$

and

$U^{2l_{2}}(p)=\displaystyle-2\,{\frac{\alpha_{2}a}{\beta_{2}}}\sqrt{\frac{1}{p}}+2\,{\frac{{a}^{2}{\alpha_{2}}^{2}}{{\beta_{2}}^{2}p}}+\,{\frac{2}{\beta_{2}}\left(\frac{1}{2}\,\alpha_{2}ab-{\frac{{\alpha_{2}}^{3}{a}^{3}}{{\beta_{2}}^{2}}}\right)\left(\frac{1}{p}\right)^{3/2}}.$

The corresponding analytical approximation in the time domain during a short
time step $t\in[0,\Delta t]$ is:

$u^{a}(l_{2},t)=u^{1l_{2}}(t)\ast g_{2}(t)+u^{2l_{2}}(t)\ast
r(l_{2},t)+2r(l_{2},t)+O\left(\Delta t\right),$

where

$u^{1l_{2}}(t)={\mathcal{L}}^{-1}\\{U^{1l_{2}}(p)\\}=-{\frac{{a}^{2}\alpha_{2}}{{\beta_{2}}^{2}}}+{\frac{a}{\sqrt{\pi\,t}\beta_{2}}}+{\frac{a\left(2\,{a}^{2}{\alpha_{2}}^{2}-b{\beta_{2}}^{2}\right)}{{\beta_{2}}^{3}}\sqrt{{\frac{t}{\pi}}}};$

and

$u^{2l_{2}}(t)={\mathcal{L}}^{-1}\\{U^{2l_{2}}(p)\\}=-2\,{\frac{\alpha_{2}a}{\sqrt{\pi\,t}\beta_{2}}}+2\,{\frac{{a}^{2}{\alpha_{2}}^{2}}{{\beta_{2}}^{2}}}+2\,{\frac{\left(-2\,{a}^{2}{\alpha_{2}}^{2}+b{\beta_{2}}^{2}\right)\alpha_{2}a}{{\beta_{2}}^{3}}\sqrt{{\frac{t}{\pi}}}}.$

    3. 3.

Likewise, for $U_{x}(l_{1},p)$ and $u_{x}(l_{1},t)$, the results are:

$U_{x}(l_{1},p)={\frac{\sqrt{b+p}\,G_{1}\left(p\right)}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}}-2\,{\frac{\alpha_{1}\sqrt{b+p}\,R\left(l_{1},p\right)}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}};$

and

$\displaystyle
U_{x}^{a}(l_{1},p)={\beta_{1}}^{-1}\,G_{1}(p)+U_{x}^{1l_{1}}(p)G_{1}(p)+U_{x}^{2l_{1}}(p)R(l_{1},p)-2\,{\frac{\alpha_{1}}{\beta_{1}}}\,R(l_{1},p)+O\left(\frac{1}{p^{2}}\right);$

where

$U_{x}^{1l_{1}}(p)=\displaystyle{\frac{a\alpha_{1}}{{\beta_{1}}^{2}}}\sqrt{\frac{1}{p}}+{\frac{{a}^{2}{\alpha_{1}}^{2}}{{\beta_{1}}^{3}p}}+{\frac{1}{\beta_{1}}\left(-\frac{1}{2}\,{\frac{a\alpha_{1}b}{\beta_{1}}}+{\frac{{a}^{3}{\alpha_{1}}^{3}}{{\beta_{1}}^{3}}}\right)\left(\frac{1}{p}\right)^{3/2}};$

and

$U_{x}^{2l_{1}}(p)=\displaystyle-2\,{\frac{{\alpha_{1}}^{2}a}{{\beta_{1}}^{2}}}\sqrt{\frac{1}{p}}-2\,{\frac{{\alpha_{1}}^{3}{a}^{2}}{{\beta_{1}}^{3}p}}-\,{\frac{2}{\beta_{1}}\left(-\frac{1}{2}\,{\frac{{\alpha_{1}}^{2}ab}{\beta_{1}}}+{\frac{{\alpha_{1}}^{4}{a}^{3}}{{\beta_{1}}^{3}}}\right)\left(\frac{1}{p}\right)^{3/2}}.$

The related analytical approximation in the time domain during a short time
step $t\in[0,\Delta t]$ is:

$u_{x}^{a}(l_{1},t)={\beta_{1}}^{-1}\,g_{1}(t)+u_{x}^{1l_{1}}(t)\ast
g_{1}(t)+u_{x}^{2l_{1}}(t)\ast
r(l_{1},t)-2\,{\frac{\alpha_{1}}{\beta_{1}}}\,r(l_{1},t)+O\left(\Delta
t\right),$

where

$u_{x}^{1l_{1}}(t)={\mathcal{L}}^{-1}\\{U_{x}^{1l_{1}}(p)\\}={\frac{a\alpha_{1}}{\sqrt{\pi\,t}{\beta_{1}}^{2}}}+{\frac{{a}^{2}{\alpha_{1}}^{2}}{{\beta_{1}}^{3}}}+{\frac{a\alpha_{1}\left(2\,{a}^{2}{\alpha_{1}}^{2}-b{\beta_{1}}^{2}\right)}{{\beta_{1}}^{4}}\sqrt{{\frac{t}{\pi}}}};$

and

$u_{x}^{2l_{1}}(t)={\mathcal{L}}^{-1}\\{U_{x}^{2l_{1}}(p)\\}=-2\,{\frac{{\alpha_{1}}^{2}a}{\sqrt{\pi\,t}{\beta_{1}}^{2}}}-2\,{\frac{{\alpha_{1}}^{3}{a}^{2}}{{\beta_{1}}^{3}}}+2\,{\frac{\left(-2\,{a}^{2}{\alpha_{1}}^{2}+b{\beta_{1}}^{2}\right){\alpha_{1}}^{2}a}{{\beta_{1}}^{4}}\sqrt{{\frac{t}{\pi}}}}.$

    4. 4.

Finally for $U_{x}(l_{2},p)$ and $u_{x}(l_{2},t)$, calculations give:

$U_{x}(l_{2},p)={\frac{\sqrt{b+p}\,G_{2}\left(p\right)}{\beta_{2}\sqrt{b+p}+\alpha_{2}a}}-2\,{\frac{\alpha_{2}\sqrt{b+p}\,R\left(l_{2},p\right)}{\beta_{2}\sqrt{b+p}+\alpha_{2}a}};$

and

$\displaystyle
U_{x}^{a}(l_{2},p)={\beta_{2}}^{-1}\,G_{2}(p)+U_{x}^{1l_{2}}(p)G_{2}(p)+U_{x}^{2l_{2}}(p)R(l_{2},p)-2\,{\frac{\alpha_{2}}{\beta_{2}}}\,R(l_{2},p)+O\left(\frac{1}{p^{2}}\right);$

where

$U_{x}^{1l_{2}}(p)=\displaystyle-{\frac{\alpha_{2}a}{{\beta_{2}}^{2}}}\sqrt{\frac{1}{p}}+{\frac{{a}^{2}{\alpha_{2}}^{2}}{{\beta_{2}}^{3}p}}+{\frac{1}{\beta_{2}}\left(\frac{1}{2}\,{\frac{\alpha_{2}ab}{\beta_{2}}}-{\frac{{\alpha_{2}}^{3}{a}^{3}}{{\beta_{2}}^{3}}}\right)\left(\frac{1}{p}\right)^{3/2}};$

and

$U_{x}^{2l_{2}}(p)=\displaystyle
2\,{\frac{{\alpha_{2}}^{2}a}{{\beta_{2}}^{2}}}\sqrt{\frac{1}{p}}-2\,{\frac{{\alpha_{2}}^{3}{a}^{2}}{{\beta_{2}}^{3}p}}-\,{\frac{2}{\beta_{2}}\left(\frac{1}{2}\,{\frac{{\alpha_{2}}^{2}ab}{\beta_{2}}}-{\frac{{\alpha_{2}}^{4}{a}^{3}}{{\beta_{2}}^{3}}}\right)\left(\frac{1}{p}\right)^{3/2}}.$

The corresponding analytical truncation in the time domain during a short time
step $t\in[0,\Delta t]$ is written as:

$u_{x}^{a}(l_{2},t)={\beta_{2}}^{-1}\,g_{2}(t)+u_{x}^{1l_{2}}(t)\ast
g_{2}(t)+u_{x}^{2l_{2}}(t)\ast
r(l_{2},t)-2\,{\frac{\alpha_{2}}{\beta_{2}}}\,r(l_{2},t)+O\left(\Delta
t\right),$

where

$u_{x}^{1l_{2}}(t)={\mathcal{L}}^{-1}\\{U_{x}^{1l_{2}}(p)\\}=-{\frac{\alpha_{2}a}{\sqrt{\pi\,t}{\beta_{2}}^{2}}}+{\frac{{a}^{2}{\alpha_{2}}^{2}}{{\beta_{2}}^{3}}}+{\frac{\left(-2\,{a}^{2}{\alpha_{2}}^{2}+b{\beta_{2}}^{2}\right)\alpha_{2}a}{{\beta_{{2}}}^{4}}\sqrt{{\frac{t}{\pi}}}};$

and

$u_{x}^{2l_{2}}(t)={\mathcal{L}}^{-1}\\{U_{x}^{2l_{2}}(p)\\}=2\,{\frac{{\alpha_{2}}^{2}a}{\sqrt{\pi\,t}{\beta_{2}}^{2}}}-2\,{\frac{{\alpha_{2}}^{3}{a}^{2}}{{\beta_{2}}^{3}}}+2\,{\frac{{\alpha_{2}}^{2}a\left(2\,{a}^{2}{\alpha_{2}}^{2}-b{\beta_{2}}^{2}\right)}{{\beta_{2}}^{4}}\sqrt{{\frac{t}{\pi}}}}.$

  * •

Case $\beta_{1}\neq 0$ and $\beta_{2}=0$

    1. 1.

As in the first case, the system is solved using the reduced determinant
$\mathrm{det}(\mathcal{S}^{a})$. This leads to:

$\displaystyle
U(l_{1},p)=-{\frac{G_{1}\left(p\right)a}{\beta_{1}\sqrt{b+p}-a\alpha_{{1}}}}+2\,{\frac{\sqrt{b+p}\,R\left(l_{1},p\right)\beta_{1}}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}},$

and $U(l_{1},p)$ and its asymptotic expansion $U^{a}(l_{1},p)$, as well as
$u(l_{1},t)$ and its truncation expansion $u^{a}(l_{1},t)$ are the same as in
the case $\beta_{1}\beta_{2}\neq 0$, item (1).

    2. 2.

Concerning $U(l_{2},p)$ and $u(l_{2},t)$, the results are reduced to:

$U(l_{2},p)={\frac{G_{2}\left(p\right)}{\alpha_{2}}};$

and the reverse in the time domain during a short time step $t\in[0,\Delta t]$
is:

$u(l_{2},t)=\frac{g_{2}(t)}{\alpha_{2}}.$

    3. 3.

The functions $U_{x}(l_{1},p)$, $u_{x}(l_{1},t)$ and their asymptotic and
truncation expansions $U^{a}_{x}(l_{1},p)$, $u^{a}_{x}(l_{1},t)$ are the same
as in the case $\beta_{1}\beta_{2}\neq 0$, item (3), since

$U_{x}(l_{1},p)={\frac{\sqrt{b+p}\,G_{1}\left(p\right)}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}}-2\,{\frac{\alpha_{1}\sqrt{b+p}\,R\left(l_{1},p\right)}{\beta_{1}\sqrt{b+p}-a\alpha_{1}}}.$

    4. 4.

About $U_{x}(l_{2},p)$ and $u_{x}(l_{2},t)$, one has:

$U_{x}(l_{2},p)={\frac{\sqrt{b+p}\,G_{2}\left(p\right)}{\alpha_{2}a}}-2\,{\frac{\sqrt{b+p}\,R\left(l_{2},p\right)}{a}};$

and the calculations will use here the assumptions that $r(l_{2},t)$ and
$g_{2}(t)$ are once differentiable relatively to the time variable $t$ on
$t\in[0,\Delta t]$, and their respective derivative $r^{\prime}(l_{2},t)$ and
$g^{\prime}_{2}(t)$ verify the relations:
$r(l_{2},t)-r(l_{2},0)=\displaystyle\int_{0}^{t}\,r^{\prime}(l_{2},\tau)d\tau$,
and $g_{2}(t)-g_{2}(0)=\displaystyle\int_{0}^{t}\,g^{\prime}_{2}(\tau)d\tau.$
The derivatives are assumed to be of exponential order so that if
$R^{q}(l_{2},p)$ and $G^{q}_{2}(p)$ are their Laplace transforms, we have
according to the properties:

$R(l_{2},p)=\displaystyle{\mathcal{L}}\left\\{\int_{0}^{t}\,r^{\prime}(l_{2},\tau)d\tau\right\\}+\frac{r(l_{2},0)}{p}=\frac{R^{q}(l_{2},p)}{p}+\frac{r(l_{2},0)}{p}$

and

$G_{2}(p)=\displaystyle{\mathcal{L}}\left\\{\int_{0}^{t}\,g^{\prime}_{2}(\tau)d\tau\right\\}+\frac{g_{2}(0)}{p}=\frac{G^{q}_{2}(p)}{p}+\frac{g_{2}(0)}{p}.$

Then

$U_{x}(l_{2},p)={\frac{\sqrt{b+p}\,\left(G^{q}_{2}(p)+g_{2}(0)\right)}{\alpha_{2}a\,p}}-2\,{\frac{\sqrt{b+p}\,\left(R^{q}(l_{2},p)+r(l_{2},0)\right)}{a\,p}};$

and

$\displaystyle
U_{x}^{a}(l_{2},p)=U_{x}^{l_{2}}(p)\left(\frac{G^{q}_{2}(p)+g_{2}(0)}{a\alpha_{2}}-2\frac{R^{q}(l_{2},p)+r(l_{2},0)}{a}\right)+O\left(\frac{1}{p^{5/2}}\right);$

where

$U_{x}^{l_{2}}(p)=\displaystyle\sqrt{\frac{1}{p}}+\,{\frac{b}{2}}\left(\frac{1}{p}\right)^{3/2}.$

In this sub-case, a truncation expansion in the time domain during a short
time step $t\in[0,\Delta t]$ is:

$u_{x}^{a}(l_{2},t)=\frac{1}{a}\left(\frac{g_{2}(0)}{\alpha_{2}}-2r(l_{2},0)\right)u_{x}^{l_{2}}(t)+\frac{1}{a}u_{x}^{l_{2}}(t)\ast\left(\frac{g^{\prime}_{2}(t)}{\alpha_{2}}-2r^{\prime}(l_{2},t)\right)+O\left(\Delta
t\right),$

where

$u_{x}^{l_{2}}(t)={\frac{bt+1}{\sqrt{\pi\,t}}}.$

  * •

Case $\beta_{1}=0$ and $\beta_{2}\neq 0$

    1. 1.

The function $U(l_{1},p)$ and its reverse $u(l_{1},t)$ in the time domain
($t\in[0,\Delta t]$) are respectively:

$U(l_{1},p)={\frac{G_{1}\left(p\right)}{\alpha_{1}}};$

and

$u(l_{1},t)=\frac{g_{1}(t)}{\alpha_{1}}.$

    2. 2.

For $U(l_{2},p)$ and $u(l_{2},t)$, the results are identical to those obtained
in the case $\beta_{1}\beta_{2}\neq 0$, item (2) since:

$U(l_{2},p)={\frac{G_{2}\left(p\right)a}{\beta_{2}\sqrt{b+p}+\alpha_{{2}}a}}+2\,{\frac{R\left(l_{2},p\right)\sqrt{b+p}\beta_{2}}{\beta_{2}\sqrt{b+p}+\alpha_{2}a}};$

    3. 3.

For $U_{x}(l_{1},p)$ and $u_{x}(l_{1},t)$, calculations give:

$U_{x}(l_{1},p)=-{\frac{\sqrt{b+p}\,G_{1}\left(p\right)}{\alpha_{1}a}}+2\,{\frac{\sqrt{b+p}\,R\left(l_{1},p\right)}{a}};$

and the assumptions are that $r(l_{1},t)$ and $g_{1}(t)$ are once
differentiable in rapport to the time variable $t$ on $t\in[0,\Delta t]$. The
respective derivatives $r^{\prime}(l_{1},t)$ and $g^{\prime}_{1}(t)$ can be
related to the original functions by the relations:
$r(l_{1},t)-r(l_{1},0)=\displaystyle\int_{0}^{t}\,r^{\prime}(l_{1},\tau)d\tau$,
and $g_{1}(t)-g_{1}(0)=\displaystyle\int_{0}^{t}\,g^{\prime}_{1}(\tau)d\tau.$
Since the derivatives are assumed to be of exponential order, their Laplace
transforms $R^{q}(l_{1},p)$ and $G^{q}_{1}(p)$ verify the properties:

$R(l_{1},p)=\displaystyle{\mathcal{L}}\left\\{\int_{0}^{t}\,r^{\prime}(l_{1},\tau)d\tau\right\\}+\frac{r(l_{1},0)}{p}=\frac{R^{q}(l_{1},p)}{p}+\frac{r(l_{1},0)}{p}$

and

$G_{1}(p)=\displaystyle{\mathcal{L}}\left\\{\int_{0}^{t}\,g^{\prime}_{1}(\tau)d\tau\right\\}+\frac{g_{1}(0)}{p}=\frac{G^{q}_{1}(p)}{p}+\frac{g_{1}(0)}{p}.$

Then

$U_{x}(l_{1},p)=-{\frac{\sqrt{b+p}\,\left(G^{q}_{1}(p)+g_{1}(0)\right)}{\alpha_{1}a\,p}}+2\,{\frac{\sqrt{b+p}\,\left(R^{q}(l_{1},p)+r(l_{1},0)\right)}{a\,p}};$

and

$\displaystyle
U_{x}^{a}(l_{1},p)=U_{x}^{l_{1}}(p)\left(-\frac{G^{q}_{1}(p)+g_{1}(0)}{a\alpha_{1}}+2\frac{R^{q}(l_{1},p)+r(l_{1},0)}{a}\right)+O\left(\frac{1}{p^{5/2}}\right);$

where

$U_{x}^{l_{1}}(p)=\displaystyle\sqrt{\frac{1}{p}}+\,{\frac{b}{2}}\left(\frac{1}{p}\right)^{3/2}.$

A truncation expansion in the time domain during a short time step
$t\in[0,\Delta t]$ corresponds to:

$u_{x}^{a}(l_{1},t)=\frac{1}{a}\left(-\frac{g_{1}(0)}{\alpha_{1}}+2r(l_{1},0)\right)u_{x}^{l_{1}}(t)+\frac{1}{a}u_{x}^{l_{1}}(t)\ast\left(-\frac{g^{\prime}_{1}(t)}{\alpha_{1}}+2r^{\prime}(l_{1},t)\right)+O\left(\Delta
t\right),$

where

$u_{x}^{l_{1}}(t)={\frac{bt+1}{\sqrt{\pi\,t}}}.$

    4. 4.

Since:

$U_{x}(l_{2},p)={\frac{\sqrt{b+p}\,G_{2}\left(p\right)}{\beta_{2}\sqrt{b+p}+\alpha_{2}a}}-2\,{\frac{\alpha_{2}\sqrt{b+p}\,R\left(l_{2},p\right)}{\beta_{2}\sqrt{b+p}+\alpha_{2}a}};$

the results for $U_{x}(l_{2},p)$ and $u_{x}(l_{2},t)$ in this sub-case are
identical to those obtained in the case $\beta_{1}\beta_{2}\neq 0$, item (4).

  * •

Case $\beta_{1}=0$ and $\beta_{2}=0$

    1. 1.

The expression of $U(l_{1},p)$ and of its reverse $u(l_{1},t)$ are identical
to those obtained in the case $\beta_{1}=0$ and $\beta_{2}\neq 0$, item (1).

    2. 2.

The expression of $U(l_{2},p)$ and of its reverse $u(l_{2},t)$ are the same as
those obtained in the case $\beta_{1}\neq 0$ and $\beta_{2}=0$, item (2).

    3. 3.

The expansions of $U_{x}(l_{1},p)$ and of its reverse $u_{x}(l_{1},t)$ are
equal to those obtained in the case $\beta_{1}=0$ and $\beta_{2}\neq 0$, item
(3).

    4. 4.

The expansions of $U_{x}(l_{2},p)$ and of its reverse $u_{x}(l_{2},t)$ are
identical those obtained in the case $\beta_{1}\neq 0$ and $\beta_{2}=0$, item
(4).

##  6 Discussion of the results and applications

First, the results mentioned in this document can be compared and discussed on
a specific example. Some powerful PDE toolbox functions exist in Matlab
software for example, and scripts based on Gauss-Seidel and finite difference
methods are available online. Therefore, the curves obtained from the Fourier
decomposition method, when using partial sums of the infinite series solution,
and those obtained from the method of approximate analytical solutions in
short time limits developed in this paper, can both be compared to numerical
solution curves, when using toolbox functions. For this purpose, the example
titled Example 6.1. is considered from [13]. It consists to solve a problem
where the initial condition matches with the boundary conditions. The equation
is stated as:

$\frac{\partial u}{\partial t}-a^{2}\frac{\partial^{2}u}{\partial
x^{2}}=0,\quad 0<x<l,\quad t>0,$ (44)

subject to consistent initial and boundary conditions

$u(x,0)=\varphi(x),\quad\quad u(0,t)=u(0,t)=0,$ (45)

where

$\varphi(x)=\left\\{\begin{array}[]{cccc}\displaystyle\frac{x}{l}u_{0}\quad\mathrm{for}\quad
0\leq x\leq\frac{l}{2}\\\ \\\
\displaystyle\frac{l-x}{l}u_{0}\quad{for}\quad\frac{l}{2}<x\leq
l\end{array}\right.$ (46)

with $u_{0}$ being a constant. The exact solution obtained by the Fourier
decomposition method can be reported as the following infinite series:

$\displaystyle
u(x,t)=\frac{4u_{0}}{\pi^{2}}\sum_{k=1}^{+\infty}\frac{(-1)^{k+1}}{(2k-1)^{2}}\exp\left(-\frac{a^{2}(2k-1)^{2}\pi^{2}}{l^{2}}t\right)\sin\frac{(2k-1)\pi
x}{l}.$ (47)

Assigning the corresponding parameter values to the general problem (1)-(4),
that is, $l_{1}=0,\ l_{2}=l,\ b=0,\ f(x,t)=0,\ \beta_{1}=\beta_{2}=0,\
\alpha_{1}=\alpha_{2}=1,\ g_{1}(t)=g_{2}(t)=0,$ the p-domain solution (27) is
reduced to:

$\begin{array}[]{ll}U(x,p)=\displaystyle\frac{a}{2\sqrt{p}}\left[U_{x}(l,p)\exp\left(\frac{-(l-x)\sqrt{p}}{a}\right)-U_{x}(0,p)\exp\left(\frac{-x\sqrt{p}}{a}\right)\right]\\\
+\displaystyle{\mathcal{L}}\left\\{\frac{1}{2a\sqrt{\pi
t}}\int_{0}^{l}\varphi(\xi)\left[\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}t}\right)\right]d\xi\right\\}.\end{array}$
(48)

Then, for the case when $l=10$, $u_{0}=5$ and $a^{2}=0.25$ as in the Example
6.1. from [13], the function $\varphi$ is rewritten as:

$\varphi(x)=\left\\{\begin{array}[]{cccc}\displaystyle\frac{1}{2}x\quad\mathrm{for}\quad
0\leq x\leq 5\\\ \\\ \displaystyle 5-\frac{1}{2}x\quad{for}\quad 5<x\leq
10;\end{array}\right.$

while the reported series (47) becomes:

$\displaystyle
u(x,t)=\frac{20}{\pi^{2}}\sum_{k=1}^{+\infty}\frac{(-1)^{k+1}}{(2k-1)^{2}}\exp\left(-\frac{(2k-1)^{2}\pi^{2}}{400}t\right)\sin\frac{(2k-1)\pi
x}{10}.$ (49)

By using equations (37) and (39) and denoting $\exp(x)$ by ${\rm e}^{x}$, the
exact operational solution (48) can be written as:

$U(x,p)=\displaystyle\frac{1}{4\sqrt{p}}\left[U_{x}(10,p)\,\exp\left(-2(10-x)\sqrt{p}\right)-U_{x}(0,p)\,\exp\left(-2x\sqrt{p}\right)\right]+R(x,p).$
(50)

In the above expression,

$U_{x}(0,p)=-U_{x}(10,p)=\displaystyle-\frac{1}{2}\,{\frac{-{{\rm
e}^{-20\,\sqrt{p}}}+2\,{{\rm e}^{-10\,\sqrt{p}}}-1}{\left({{\rm
e}^{-20\,\sqrt{p}}}+1\right)p}};$ (51)

and

$R(x,p)=\left\\{\begin{array}[]{cccc}\displaystyle-\frac{1}{8}\,{\frac{-4\,x\sqrt{p}+2\,{{\rm
e}^{2\,\left(-5+x\right)\sqrt{p}}}-{{\rm e}^{-2\,x\sqrt{p}}}-{{\rm
e}^{2\,\left(-10+x\right)\sqrt{p}}}}{{p}^{\frac{3}{2}}}}\quad\mathrm{for}\quad
0\leq x\leq 5,\\\ \\\
\displaystyle-\frac{1}{8}\,{\frac{4\,x\sqrt{p}-40\,\sqrt{p}+2\,{{\rm
e}^{-2\,\left(-5+x\right)\sqrt{p}}}-{{\rm e}^{-2\,x\sqrt{p}}}-{{\rm
e}^{2\,\left(-10+x\right)\sqrt{p}}}}{{p}^{\frac{3}{2}}}}\quad{\rm for}\quad
5<x\leq 10.\end{array}\right.$

The formula (23) can then be checked, and the exactness of the operational
solution (50) is proved, that is:

$-a^{2}\frac{d^{2}U}{dx^{2}}(x,p)+pU(x,p)=\varphi(x).$

Now, according to the results obtained in section 5 relatively to the present
case ($\beta_{1}=\beta_{2}=0$), the corresponding truncation expansion in the
time domain during the short time step $t\in[0,\Delta t]$, are respectively
recalled as:

$u^{a}(x,t)=r(x,t)+o\left((\Delta t)^{n}\right);$

and

$u_{x}^{a}(l_{1},t)=\frac{1}{a}\left(-\frac{g_{1}(0)}{\alpha_{1}}+2r(l_{1},0)\right)u_{x}^{l_{1}}(t)+\frac{1}{a}u_{x}^{l_{1}}(t)\ast\left(-\frac{g^{\prime}_{1}(t)}{\alpha_{1}}+2r^{\prime}(l_{1},t)\right)+O\left(\Delta
t\right),$

where

$u_{x}^{l_{1}}(t)={\frac{bt+1}{\sqrt{\pi\,t}}},$

and again

$u_{x}^{a}(l_{2},t)=\frac{1}{a}\left(\frac{g_{2}(0)}{\alpha_{2}}-2r(l_{2},0)\right)u_{x}^{l_{2}}(t)+\frac{1}{a}u_{x}^{l_{2}}(t)\ast\left(\frac{g^{\prime}_{2}(t)}{\alpha_{2}}-2r^{\prime}(l_{2},t)\right)+O\left(\Delta
t\right),$

where

$u_{x}^{l_{2}}(t)={\frac{bt+1}{\sqrt{\pi\,t}}}.$

By using Maple software for example, the analytical approximations in short
time limits ($t\in[0,\Delta t]$), can be computed from above expressions as:

$\begin{array}[]{cccc}u^{a}(x,t)\simeq r(x,t)=\displaystyle\frac{1}{\sqrt{\pi
t}}\int_{0}^{10}\varphi(\xi)\left[\exp\left(-\frac{(\xi-x)^{2}}{t}\right)\right]d\xi\\\
\displaystyle=\frac{1}{4}\left((x-10){\rm
erf}\left(\frac{-10+x}{\sqrt{t}}\right)+(-2x+10){\rm
erf}\left(\frac{-5+x}{\sqrt{t}}\right)+x\>{\rm
erf}\left(\frac{x}{\sqrt{t}}\right)\right)\\\
\displaystyle+\frac{\sqrt{t}}{4\sqrt{\pi}}\left(-2\,{{\rm
e}^{-{\frac{\left(x-5\right)^{2}}{t}}}}+{{\rm
e}^{-{\frac{\left(x-10\right)^{2}}{t}}}}+{{\rm
e}^{-{\frac{{x}^{2}}{t}}}}\right);\end{array}$ (52)

and

$u_{x}^{a}(0,t)=-u_{x}^{a}(10,t)=-\frac{1}{2}{\rm
erf}\left(\frac{10}{\sqrt{t}}\right)+{\rm
erf}\left(\frac{5}{\sqrt{t}}\right)+O\left(\Delta t\right)$ (53)

where ${\rm erf}(x)$ is the Error Function defined by:

$\displaystyle{\rm
erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\mathrm{e}^{y^{2}}dy.$

Figure 1: (a) Solution $u(x,t)$ at $t=\Delta t=10^{-2}.$ (b) Partial
derivative $u_{x}(x,t)$ at $t=\Delta t=10^{-2}$

Compared to the function $\varphi(x)$ characterizing the initial condition,
it’s remarkable that the derivatives coincide at $l_{1}=0$ as well as
$l_{2}=10$ , that is :

$\displaystyle\lim_{t\rightarrow
0}u_{x}^{a}(0,t)=\varphi^{\prime}(0)=\frac{1}{2},$ (54)

and

$\displaystyle\lim_{t\rightarrow
10}u_{x}^{a}(10,t)=\varphi^{\prime}(10)=-\frac{1}{2}.$ (55)

In the limits of our knowledge, the truncation expansion (53) was not
available directly from the p-domain formula (51).

Figure 1 shows curves of the solution of the problem (44)-(46) and of its
derivative during the first time step $t\in[0,\Delta t]$ with $\Delta
t=10^{-2}$. The curves are obtained by using the three different methods. The
numerical solution is computed by using toolbox functions of Matlab, while the
series solution provided by the Fourier decomposition method is truncated at
its first 20 terms. The approximate analytical solution is represented by the
function $u^{a}$ given by the above formula (52). The three resulting curves
of the solution $u$ are almost identical as shown on figure 1(a). Their
shapes, except around the peak point of abscissa $x=5$, seem to be very
similar to that of the initial condition function (46), since the time passed
from $t=0$ to $t=\Delta t=10^{-2}$ is still relatively small. On the figure
1(b), except at the domain boundaries, i.e. at $x=0$ and $x=10$, the
derivative curves almost coincide for the approximate analytical and the
numerical solutions. However, the curve for the series solution exhibits
deviations in the form of small oscillations along the two former curves,
especially around the peak point ($x=5$). This may suggest on the one hand,
the potential consistency of the approximate analytical solution as well as
that of the numerical method to properly accounting for nonphysical phenomena
as derivative jumps and the so-called Gibbs phenomenon arising in generalized
series solutions. On the other hand, this may illustrate the pertinence of
refined truncation expansions at the domain boundaries for the solution, that
is for $u^{a}(0,t)$ and $u^{a}(10,t)$ which are keep null here, and for its
derivative, i.e., $u^{a}_{x}(0,t)$ and $u^{a}_{x}(10,t)$ in the present case.
Figures 2(a) and 2(b) respectively show evolution of the derivatives
$u^{a}_{x}(0,t)$ and $u^{a}_{x}(10,t)$ from formulae (53). During the short
time step $[0,\Delta t]$, the time-varying derivative curves obtained from the
approximate and series solutions can be compared. While the convergence of the
series solution is poor especially at points close to the ends of the domain,
the shapes shown by the curves of $u^{a}_{x}(0,t)$ and $u^{a}_{x}(10,t)$
correspond well to the similitude highlighted above with the initial function
and confirmed by the calculations of the limits (54) and (55). Beyond its
precision or consistency, another advantage of performing the approximate
analytical solutions $u^{a}_{x}(0,t)$ and $u^{a}_{x}(10,t)$ during the short
time limit $[0,\Delta t]$, consists of the computational efficiency of these
solutions. Indeed, at the first execution of the code source, the script
running time of their procedure is far reduced compared to that of the series
solution (about 20 times).

Figure 2: (a) Evolution of the derivative $u_{x}(0,t)$ for $t\in[0,\Delta t]$,
$\Delta t=10^{-2}$. (b) Evolution of the derivative $u_{x}(10,t)$ for
$t\in[0,\Delta t]$, $\Delta t=10^{-2}$

Next, the exact solution in the Laplace domain can be extended to unbounded
domains as infinite or semi-infinite intervals, for the space variable $x$.
Indeed when $l_{1}=-\infty$ and $l_{2}=+\infty$, one has
$U_{x}(l_{1},p)=U_{x}(l_{2},p)=0$, and by taking the limits in equation (27),
the exact p-domain solution is reduced to $U(x,p)=R(x,p)$, which corresponds
to the following solution in the time domain:

$\begin{array}[]{ll}u(x,t)=r(x,t)=\displaystyle\frac{\exp(-bt)}{2a\sqrt{\pi
t}}\int_{-\infty}^{+\infty}\varphi(\xi)\left[\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}t}\right)\right]d\xi\\\
+\displaystyle\frac{1}{2a\sqrt{\pi}}\int_{0}^{t}d\theta\int_{-\infty}^{+\infty}\frac{\exp(-b(t-\theta))}{\sqrt{(t-\theta)}}\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}(t-\theta)}\right)f(\xi,\theta)d\xi.\end{array}$

In the case when $b=0$, the above solution is identical to the one reported by
[13] in their book at the section 6.8 titled: "The Heat Equation in an
Infinite Region", where the Fourier decomposition method was used. About semi-
infinite domains, let us chose $l_{1}=0$ and $l_{2}=+\infty$ for example, then
$U_{x}(l_{2},p)=0$, and $u(l_{2},t)$ is to be considered as constant. The
p-domain solution (27) reduces to:

$U(x,p)=\displaystyle\frac{1}{2}U(0,p)\exp\left(\frac{-x\sqrt{b+p}}{a}\right)-\frac{aU_{x}(0,p)}{2\sqrt{b+p}}\exp\left(\frac{-x\sqrt{b+p}}{a}\right)+R(x,p),$
(56)

where $R(x,p)={\mathcal{L}}\\{r(x,t)\\}$, and

$\begin{array}[]{ll}r(x,t)=\displaystyle\frac{\exp(-bt)}{2a\sqrt{\pi
t}}\int_{0}^{+\infty}\varphi(\xi)\left[\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}t}\right)\right]d\xi\\\
+\displaystyle\frac{1}{2a\sqrt{\pi}}\int_{0}^{t}d\theta\int_{0}^{+\infty}\frac{\exp(-b(t-\theta))}{\sqrt{(t-\theta)}}\exp\left(-\frac{(\xi-x)^{2}}{4a^{2}(t-\theta)}\right)f(\xi,\theta)d\xi.\end{array}$
(57)

Setting as in the statement of the Problem a. of the chapter 8.1 in [20],
$b=0$, $\varphi(x)=t_{0}=\mathrm{(constant)}$,
$u(0,t)=g_{1}(t)=t_{a}=\mathrm{(constant)}$ implying $U(0,p)=t_{a}/p$, it can
be verified for $f(x,t)={w}/{c\gamma}=\mathrm{(constant)}$ that:

$Ux(0,p)=\displaystyle\frac{2\sqrt{p}}{a}R(0,p)-\frac{t_{a}}{a\sqrt{p}};$

and,

$R\left(x,p\right)=\displaystyle{\frac{t_{0}}{p}}+{\frac{w}{{p}^{2}c\gamma}}-\,{\frac{t_{0}}{2p}{\exp\left({-{\frac{x\sqrt{p}}{a}}}\right)}}-\,{\frac{w}{2{p}^{2}c\gamma}{\exp\left({-{\frac{x\sqrt{p}}{a}}}\right)}}.$

This leads to reduce the exact p-domain solution (56) into the form of:

$U\left(x,p\right)=\displaystyle{\frac{t_{0}}{p}}+{\frac{w}{{p}^{2}c\gamma}}+\,{\frac{(t_{a}-t_{0})}{p}{\exp\left({-{\frac{x\sqrt{p}}{a}}}\right)}}-\,{\frac{w}{{p}^{2}c\gamma}{\exp\left({-{\frac{x\sqrt{p}}{a}}}\right)}};$

which is identical to the solution (8.1.11) reported in [20], provided that
$a$ is replaced by $\sqrt{a}$ as specified in the statement of the problem.
Note in the latter reference that the Laplace transform method was directly
used to solve the problem, since the initial condition is specified as a
constant function.

Finally, as a perspective for numerical simulation purposes, schemes to be
used for deriving the complete analytical approximation of the solution
$u(x,t)$ on $[l_{1},l_{2}]\times[0,T]$ may be obtained. The procedure for
calculating analytical approximations during the first time step
$t\in[0,\Delta t]$ can be repeated for any short time step by updating the
source term $f(x,t)$, the functions $\varphi(x)$, $g_{1}(t)$ and $g_{2}(t)$
that are related to the initial and boundary conditions. Also, the largest
class to which the above functions must belong in order to obtain the exact
solution in the Laplace domain may be studied more precisely. Moreover,
computational models for one-dimensional Stefan problems as reviewed in [16]
can be efficiently handled when using the truncation expansions
$u^{a}(l_{1},t)$, $u^{a}_{x}(l_{1},t)$, $u^{a}(l_{2},t)$ and
$u^{a}_{x}(l_{2},t)$ of the approximate analytical solution at the boundaries
of the domain as expressed in section 5. An example of using alike formulae
for a specific problem of a spherically symmetric droplet evaporation can be
seen in [1].

##  7 Conclusion

This study has permitted to calculate an explicit solution in the Laplace
domain and analytical approximations in the earlier time-step to initial
boundary value problems for the one-dimensional parabolic equation with
constant coefficients. The problem is solved in its most general form, with
the boundary conditions stated in a unified way on any bounded generic
interval of the real line. Compared to the classical or generalized series
solutions that can be obtained for the same problem by using the Fourier
decomposition method, the analytical approximations in short time limits are
proven to be more consistent and sufficiently simple to improve computational
efficiency in numerical schemes and simulations. Early time behaviors of heat
or mass reaction-diffusion processes are of great interest in engineering and
have a wide range of applications in fields like Computational Fluid Dynamics
(CFD) and Nuclear Energy. In addition, the exact operational solution obtained
for the problem can be extended to unbounded domains. This explicit solution
can also be considered as an actual advance in the study of linear parabolic
equations. Though they may be sought by means of Laplace inversion theorems,
the number of exact inverse transforms from the Laplace domain into the time
domain are limited in most tables or by using software facilities. However,
inverse Laplace transforms can be accomplished numerically regardless of the
complexity of the operational solutions. Therefore, accurate curves of the
complete solution in the time domain can be obtained from the exact
operational solution, when using a numerical inverse Laplace transform.

## References

  * [1] Anani, K. (2021). An efficient approximate analytical model for droplets transient heating and evaporation, Int. J Numer. Methods Appl., Vol. 20, No. 2, pp. 157-172.
  * [2] Anani, K. (2022). Series solutions for the spherically symmetric droplet transient heating problem, AIP Conference Proceedings, Vol. 2668, No. 1, pp. 050001-1–050001-8.
  * [3] Bougoffa, L., Rach, R., Wazwaz, A. M., and Duan, J. S. (2015). On the Adomian decomposition method for solving the Stefan problem, International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 25, No. 4, pp. 912–928.
  * [4] Brenn, G. (2017). _Analytical Solutions for Transport Processes. Mathematical Engineering_ , Springer, Berlin.
  * [5] Carr, E.J. and March, N.G. (2018). Semi-analytical solution of multilayer diffusion problems with time-varying boundary conditions and general interface conditions, Appl. Math. Comput., Vol. 333, pp. 286–303.
  * [6] Chang, J.Y., Chen, R.Y. and Tsai, C.C. (2022). Hermite Method of Approximate Particular Solutions for Solving Time-Dependent Convection-Diffusion-Reaction Problems, Mathematics, Vol. 10(2), No. 188, pp. 1–20.
  * [7] Cioranescu, D., Donato, P. and Roque, M.P. (2018). _An Introduction to Second Order Partial Differential Equations: Classical and Variational Solutions_ , World Scientific, Singapore.
  * [8] Debnath, L. and Bhatta, D. (2014). _Integral Transforms and Their Applications_ , CRC press, Boca Raton.
  * [9] Dobrushkin V. (2017). _Applied Differential Equations with Boundary Value Problems_ , CRC Press, Boca Raton.
  * [10] Ghasemia, M. and Kajani, M.T. (2010). Applications of He’s homotopy perturbation method to solve a diffusion convection problem, Math. Sci., Vol. 4 No. 2, pp. 171–186.
  * [11] Han J.C. (2016). _Analytical Heat Transfer_ , CRC Press, Boca Raton.
  * [12] Henner, V., Belozerova, T. and Forinash, K. (2009). _Mathematical Methods in Physics: Partial Differential Equations, Fourier Series, and Special Functions_ , CRC Press, Boca Raton.
  * [13] Henner, V., Belozerov, T. and Nepomnyashchy, A. (2019). _Partial Differential Equations: Analytical Methods and Applications_ , CRC Press, Boca Raton.
  * [14] Herron I.H. and Foster, M.R. (2008). _Partial Differential Equations in Fluiddynamics_ , Cambridge University Press, Cambridge.
  * [15] Izadi, M. and Yuzbasi, S. (2022). A hybrid approximation scheme for 1-D singularly perturbed parabolic convection-diffusion problems, Mathematical Communications, Vol. 27, No. 1, pp. 47–62.
  * [16] Javierre, E., Vuik, C., Vermolen, F. J. and Van der Zwaag, S. (2006). A comparison of numerical models for one-dimensional Stefan problems, Journal of Computational and Applied Mathematics, Vol. 192 No. 2, pp. 445–459.
  * [17] Li, X.Y., Wang, H.L. and Wu, B.Y. (2022). A stable and efficient technique for linear boundary value problems by applying kernel functions. Applied Numerical Mathematics, Vol. 172, pp. 206–214.
  * [18] Lin, J., Xu, Y. and Zhang, Y. (2020). Simulation of linear and nonlinear advection–diffusion–reaction problems by a novel localized scheme. Applied Mathematics Letters, Vol. 99, 106005.
  * [19] Liu, Y. and Zhao, X. (2010). He’s variational iteration method for solving Convection diffusion equations, Adv. Intell. Comput. Theories Appl., Vol. 6215, pp. 246–251.
  * [20] Luikov A.V. (2012). _Analytical Heat Diffusion Theory_ , Academic Press, New York.
  * [21] Marin, M., and Öchsner, A. (2018). _Essentials of Partial Differential Equations_ , Springer International Publishing, Switzerland.
  * [22] Minkowycz, W.J., Sparrow, E.M. and Murthy, J.Y. (2006). _Handbook of Numerical Heat Transfer_ , John Wiley & Sons, New York.
  * [23] Momani, S. (2007). An algorithm for solving the fractional convection diffusion equation with nonlinear source term. Commun. Nonlinear Sci. Numer. Simul., Vol. 12, No. 7, pp. 1283–1290.
  * [24] Poularikas, A.D. (2018). _Laplace Transforms, the Handbook of Formulas and Tables for Signal Processing_ , CRC press, Boca Raton.
  * [25] Su, L.D., Jiang, Z W. and Jiang, T.S. (2015). Numerical method for one-dimensional convection-diffusion equation using radical basis functions, J. Phys. Math., Vol. 6, No. 1, pp. 1–5.
  * [26] Tassaddiq, A., Yaseen, M., Yousaf, A. and Srivastava, R. (2021). Computational study of the convection-diffusion equation using new cubic B-spline approximations, AIMS Mathematics, Vol. 6, No. 5, pp. 4370–4393.
  * [27] Urban, K. and Patera, A. (2014). An improved error bound for reduced basis approximation of linear parabolic problems. Mathematics of Computation, Vol. 83, No. 288, pp. 1599–1615.
  * [28] Yüzbaşı, Ş. and Şahin, N. (2013). Numerical solutions of singularly perturbed one-dimensional parabolic convection-diffusion problems by the Bessel collocation method, Appl. Math. Comput., Vol. 220 pp. 305–315.
  * [29] Zheng, L. and Zhang, X. (2017). _Modeling and Analysis of Modern Fluid Problems_ , Academic Press, New York.