Skip to content
Publicly Available Published by De Gruyter July 31, 2020

Initial-boundary value problems for the one-dimensional linear advection–dispersion equation with decay

  • Guenbo Hwang EMAIL logo

Abstract

Initial-boundary value problems for the one-dimensional linear advection–dispersion equation with decay (LAD) are studied by utilizing a unified method, known as the Fokas method. The method takes advantage of the spectral analysis of both parts of Lax pair and the global algebraic relation coupling all initial and boundary values. We present the explicit analytical solution of the LAD equation posed on the half line and a finite interval with general initial and boundary conditions. In addition, for the case of periodic boundary conditions, we show that the solution of the LAD equation is asymptotically t-periodic for large t if the Dirichlet boundary datum is periodic in t. Furthermore, it can be shown that if the Dirichlet boundary value is asymptotically periodic for large t, then so is the unknown Neumann boundary value, which is uniquely characterized in terms of the given asymptotically periodic Dirichlet boundary datum. The analytical predictions for large t are compared with numerical results showing the excellent agreement.

1 Introduction

In this paper, we study initial-boundary value problems for the one-dimensional linear advection–dispersion equation with decay (LAD) [1]

(1)qt=Dqxxcqxrcq,

where D > 0, c ≥ 0 and rc ≥ 0 represents the dispersion, advection coefficients and the first order decay rate. This equation is also known as the convection–diffusion equation with a reaction, according to scientific contexts [2]. The LAD equation is a combination of the linear dispersion and advection equations plus the decay, but which is also an analytical model, describing physical, chemical or biological diffusive transport phenomena in science and engineering [3], [4], [5], [6]. For example, it describes a mathematical model for one-dimensional contaminant transport in porous medium systems with first-order decay. The LAD equation without the decay term was formally solved in [7], [8]. Analytical solutions of the LAD equation, including the decay term, have been developed in several literatures [1], [3], [6] (see also [9] and references therein). It should be noted that most proposed analytical solutions are for the LAD equation with relatively limited initial and boundary conditions, such as homogeneous, constant or exponential boundary conditions in time. Moreover, these solutions can be involved with the complementary error functions. A few explicit analytical solutions are available if more general or complicated initial and boundary conditions are prescribed.

We propose here explicit analytical solutions for initial-boundary value problems of the LAD equation with general boundary conditions by utilizing a unified method, also known as the Fokas method [10], [11] (see the monograph [12] and the pedagogical literature of the method [13]). The Fokas method has been introduced to analyze boundary value problems for nonlinear integrable systems, which can be considered as a significant extension of the inverse scattering transform. Moreover, it has been shown that the method can be extensively applicable to a large class of partial differential equations (PDEs); for example, nonlinear integrable systems and linear evolution equations [11], [14], [15], [16] and linear and nonlinear elliptic PDEs [17], [18], [19], [20], [21] as well as linear and integrable nonlinear discrete equations [22], [23].

The presented method has several advantages. For linear cases, the Fokas method is relatively simple, but very effective, to implement. The main steps of the Fokas method can be summarized as follows: (i) simultaneous analysis of the both parts of Lax pair, which can be considered a novel type of separability [24]; (ii) analysis of the global relation which is an algebraic equation that involves all initial and boundary values. This global relation can be used to determine the unknown boundary values, known as the generalized Dirichlet-to-Neumann map [25] (see also [26], [27] for the recent application of the global relation). Moreover, the Fokas method presents an integral representation of the solution with explicit exponential (x, t)-dependence. Thus, it enables to characterize the asymptotic behaviors of the solution and the unknown boundary values for large t [28], [29], [30], [31], [32], [33].

The purpose of the work is to present new explicit analytical solutions for initial-boundary value problems of the LAD equation with general initial and boundary conditions. More specifically, letting t=t/D and dropping prime for simplicity, we consider the following rescaled equation

(2)qt=qxxbqxaq,a,b0.

Note that one can reduce LAD equation (2) to the LAD equation without decay term by changing variable q(x,t)=eatu(x,t). The equation without decay term has been studied and solved by using the Fokas method in [34]. However, the form of equation (2) has the advantage of studying periodic asymptotics of the solution and boundary values for large t. We will derive the explicit integral representation of the solution for the Dirichlet boundary value problem of equation (2) posed on half line via the Fokas method. In addition, we will discuss the asymptotic behaviors of the solution and boundary values for large t. We will show that if the Dirichlet boundary value is periodic, the solution of the LAD is asymptotically t-periodic as t → . Moreover, it can be shown that if the Dirichlet boundary value is asymptotically periodic as t → , then so is the unknown Neumann boundary value, which can be uniquely characterized in terms of the asymptotically periodic Dirichlet boundary value. It should be noted that the Fokas method can be effectively applicable to solve boundary value problems with more general and complicated boundary conditions. In this respect, we will further address the explicit solutions for equation (2) posed on the half line with the Neumann and Robin boundary conditions as well as the solution of equation (2) formulated on a finite interval. In addition, we will present several examples including comparison of the analytical and numerical results, as applications.

2 The LAD equation on the half line

In Section 2, we first study the Dirichlet boundary value problem for the LAD equation posed on the half line

(3)qt=qxxbqxaq,x>0,t>0,

with the initial condition q(x,0)=q0(x) and the Dirichlet boundary condition, denoted by

q(0,t)=g0(t).

Also, we denote the unknown Neumann boundary value by qx(0,t)=g1(t). We assume that q0(x) is sufficiently smooth and decays rapidly as x →  and g0(t) is sufficiently smooth.

The LAD equation can be written as an overdetermined linear system, known as a Lax pair

(4a)μxikμ=q,
(4b)μt+ω(k)μ=(ikb)q+qx,

where μ=μ(x,t,k) is the eigenfunction with the spectral parameter k and the dispersion relation is given by

ω(k)=k2+ibk+a.

Indeed, the LAD equation (2) is the compatible condition of the Lax pair equation (4) in the sense that μtx=μxt implies that q solves equation (2) if the spectral parameter k is independent of x and t. Note that Reω(k)=k12k22bk2+a for k=k1+ik2 (k1,k2). We introduce the region D={k|Reω(k)<0} and we decompose the region D into D=D+D, where D+={k|Reω(k)<0 and Imz>0} and D={k|Reω(k)<0 and Imz<0} (see Figure 1). The Lax pair equation (4) can be written in the divergence form

(5)(eikx+ω(k)tq)t=(eikx+ω(k)tQ)x,

where Q(x,t,k)=(ikb)q+qx. To analyze equation (5), we introduce the Fourier transform pair:

q^(k,t)=0dxeikxq(x,t),Imk0,

and

q(x,t)=12πdkeikxq^(k,t).
Figure 1: The region D=D+∪D−$D={D}\hat{+}\cup {D}\hat{-}$ (shaded) in the complex k-plane, where Reω(k) < 0 with ω(k)=k2+ibk+a$\omega \left(k\right)={k}\hat{2}+ibk+a$ and the vertices at k=i(−b±b2+4a)/2$k=i\left(-b{\pm}\sqrt{{b}\hat{2}+4a}\right)/2$. The steepest descent contour Cs${C}_{s}$ passing through the point k = −ib/2. Efficient contour L=L1∪L2$L={L}_{1}\cup {L}_{2}$ for numerical scheme (see the text for some details).
Figure 1:

The region D=D+D (shaded) in the complex k-plane, where Reω(k) < 0 with ω(k)=k2+ibk+a and the vertices at k=i(b±b2+4a)/2. The steepest descent contour Cs passing through the point k = −ib/2. Efficient contour L=L1L2 for numerical scheme (see the text for some details).

Then, taking 0dx in equation (5) yields

(6)(eω(k)tq^(k,t))t=eω(k)tQ(0,t,k).

Let q^(k,0)=q^0(k) and we define the t-transform as

(7)g^j(k,t)=0tdseω(k)sxjq(0,s),j=0,1.

Applying 0tds in equation (6), we find the global relation given by

(8)eω(k)tq^(k,t)=q^0(k)Q^(k,t),Imk0,

where

(9)Q^(k,t)=(ikb)g^0(k,t)+g^1(k,t).

Employing the inverse Fourier transform in equation (8), the solution q(x, t) can be recovered as

(10)q(x,t)=12πdkeikxω(k)tq^0(k)12πdkeikxω(k)tQ^(k,t).

Since the integrand in the second integral of equation (10) is analytic and bounded in the region, where Im > 0 and Reω(k) > 0, by the Cauchy theorem, we deform the contour (−, ) to D+ (see Figure 1). Thus, we find the reconstruction formula for the solution q(x, t) as

(11)q(x,t)=12πdkeikxω(k)tq^0(k)12πD+dkeikxω(k)tQ^(k,t).

Note that the representation of the solution equation (11) involves the unknown function g^1(k,t) via the unknown Neumann boundary value qx(0, t). Thus, for the explicit representation of the solution, it is necessary to determine this unknown boundary value. To this end, we use the global relation equation (8), which couples all initial and boundary values.

Note that if ω(k)=ω(ν(k)), the transformation kν(k) leaves ω(k) invariant. As a consequence, g^j(k,t) (j=0,1) is invariant under this transformation. The equation ω(k)=ω(ν(k)) has a trivial root ν(k)=k and a nontrivial root ν(k)=kib. Note that if kD+, then kibD. Hence, replacing kkib in equation (8) and solving the resulting equation for g^1(k,t), the unknown boundary value is given by

(12)g^1(k,t)=q^0(kib)+ikg^0(k,t)eω(k)tq^(kib,t).

Substituting equation (12) into equation (9), we find

(13)Q^(k,t)=q^0(kib)+(2ikb)g^0(k,t)eω(k)tq^(kib,t).

The right-hand-side of equation (13) contains the unknown function q^(kib,t), which is the transform of the solution q(x, t). However, when this term is inserted into equation (11), the integrand is analytic and bounded in D+. Hence, the integral of the term involving q^(kib,t) over D+ vanishes by the Cauchy theorem. Finally, the explicit representation for the solution of equation (3) is given by

(14)q(x,t)=12πdkeikxω(k)tq^0(k)12πD+dkeikxω(k)t[q^0(kib)+(2ikb)g^0(k,t)].

Remark 2.1.

We can derive the unknown Neumann boundary value directly from equation (12) as

(15)g1(t)=12iπdkkeω(k)tq^0(k)+12iπD+dkkeω(k)t[q^0(kib)+(2ikb)g^0(k,t)].

Indeed, the effective portion of g^1(k,t) in equation (12) is given by g^1(k,t)=q^0(kib)+ikg^0(k,t). Taking the inversion formula for the t-transform equation (a.3) discussed in Appendix, we find

(16)g1(t)=12iπD+dkeω(k)t(2k+ib)q^0(kib)+12iπD+dkeω(k)tk(2ikb)q^0(kib).

Note that letting kkib (and then kibD for kD+), we find

12iπD+dkeω(k)t(k+ib)q^0(kib)=12iπDdkkeω(k)tq^0(k).

Then, we deform the contour D to (−, ) in the negative direction by using the Cauchy theorem. Thus, substituting the resulting equation into the first integral of equation (16), we find the unknown Neumann value g1(t) given by equation (15). Note that equation (15) is identical with the expression obtained from differentiating equation (14) with respect to x and evaluating the resulting equation at x = 0.

The representation of the solution given in equation (14) has the explicit (x, t)-dependence of the exponential form. Thus, it is possible to study the appropriate asymptotics of the solution for large t. Below, we will show that the solution for the LAD equation with a periodic Dirichlet boundary datum is asymptotically periodic for large t. The similar result of the Proposition 2.1 also can be found in [33] for the linearized Korteweg–de Vries equation.

Proposition 2.1.

Let q(x, t) be the solution given by equation (14). Assume that q0(x)S[0,) and g0(t)L[0,), where S denote the space of Schwartz functions. Suppose that g0(t) is a periodic function with period τ. Then for any x > 0,

(17)limnq(x,nτ+t)=q(x,t),

where q(x,t) is a periodic function with period τ and q(0,t)=g0(t). In particular, for any x > 0, q(x, t) is asymptotically periodic such that

(18)q(x,τ+t)q(x,t)=O(eat),t.

Proof.

Note that eω(k)tq^0(kib) is analytic and bounded in the region, where Im> 0 and Reω(k) > 0. Also note that by integration by parts, (2ikb)g^0(k,t)=O(1/k) as k → . Thus by the Cauchy theorem, we can deform the contour D+ to (−, ) in the second integral of equation (14) and hence, we write the representation of the solution given in equation (14) as the sum of three integrals q(x,t)=Q(1)(x,t)+Q(2)(x,t)+Q(3)(x,t), where

Q(1)(x,t)=12πdkeikxω(k)tq^0(k),Q(2)(x,t)=12πdkeikxω(k)tq^0(kib),Q(3)(x,t)=12πdkeikxω(k)t(2ikb)g^0(k).

We will estimate each integral in the above equations. For Q(1)(x,t), we write

Q(1)(x,t)=12π0dyq0(y)dkeik(xy)ω(k)t.

Then, using the following identity,

(19)dkeikxω(k)t=πtexp[at(xbt)24t],

we find

(20)|Q(1)(x,t)|M1eattq0L1[0,)

for some constant M1 > 0. For Q(2)(x,t), we write

Q(2)(x,t)=12π0dyq0(y)dkeikx+(ikb)yω(k)t.

From the identity,

dkeikx+(ikb)yω(k)t=πtexp[at(xbt)2+2(x+bt)y+y24t].

it follows that

(21)|Q(2)(x,t)|M2eattq0L1[0,)

for some constant M2 > 0.

Regarding Q(3)(x,t), we write Q(3)(x,t) as

Q(3)(x,t)=12π0tdsg0(s)Φ(x,ts),

where

Φ(x,t)=dkeikxω(k)t(2ikb).

Note that since ω(k)=2k+ib,

Φ(x,t)=1itdkeikxddk(eω(k)t)=xtdkeikxω(k)t.

Then using equation (19), we find

(22)Φ(x,t)=xπt3/2exp[at(xbt)24t].

In order to derive equation (17), let {qn(x,t)}1 be the sequence given by

(23)qn(x,t)=q(x,nτ+t)=Qn(1)(x,t)+Qn(2)(x,t)+Qn(3)(x,t),

where Qn(j)(x,t)=Q(j)(x,nτ+t) for j = 1, 2, 3. From equations (20) and (21), it follows that

|Qn(1)(x,t)+Qn(2)(x,t)|(M1+M2)ea(nτ+t)nτ+tq0L1[0,).

Hence, we find

limn[Qn(1)(x,t)+Qn(2)(x,t)]=0.

For Qn(3)(x,t), we consider the following series

(24)m=1[Qm+1(3)(x,t)Qm(3)(x,t)].

Since g0(t) is a periodic function with period τ, we know that

Qm+1(3)(x,t)Qm(3)(x,t)=12πmτ+t(m+1)τ+tdsg0(ts)Φ(x,s),

which implies that by equation (22),

|Qm+1(3)(x,t)Qm(3)(x,t)|M3ea(mτ+t)(mτ+t)3/2g0L[0,)

for some constant M3 > 0. Thus, the series given in equation (24) converges and so does the sequence Qn(3)(x,t) as n. Thus, the limit of the sequence {qn(x,t)} given in equation (23) exists for any x > 0 and we denote this limit by q(x,t). Note that

q(x,τ+t)q(x,t)=[q(x,τ+t)qn(x,τ+t)]+[qn+1(x,t)q(x,t)].

Then, taking n →  in the right-hand-side of the above equation, we find q(x,τ+t)=q(x,t). Thus, q(x,t) is a periodic function with period τ. Moreover, qn(0,t)=g0(t) yields q(0,t)=g0(t). □

Next, we study the asymptotically periodic boundary data for large t. We assume that the Dirichlet boundary value is asymptotically periodic for large t in the sense that [29]

g0(t)=g˜0(t)+O(t7/2),t,

where g˜0(t) is a periodic function with period τ. Without loss of generality, we consider the case of τ = 2π. In what follows, we will show that the unknown Neumann boundary value qx(0,t) is also asymptotically periodic for large t with the same period τ.

We first note that by integration by parts, we write the function g^0(k,t) as

g^0(k,t)=1ω(k)(eω(k)tg0(t)g0(0))1ω(k)0tdseω(k)sg˙0(s)ds.

Deforming the contour D+ to pass to the right of the singularity k=i2(b+b2+4a), which is denoted by D0+, Equation (14) can be written as

(25)q(x,t)=12πdkeikxω(k)tq^0(k)12πD+dkeikxω(k)tq^0(kib)+g0(t)+12πD0+dkeikxω(k)t2ikbω(k)×[g0(0)+0tdseω(k)sg˙0(s)],

where we have used the identity: by the residue theorem,

12πD0+dkeikx2ikbω(k)g0(t)=g0(t).

Then, differentiating equation (25) with respect to x and evaluating the resulting equation at x = 0, the function g1(t) is given by

(26)g1(t)=12iπdkkeω(k)tq^0(k)+12iπD+dkkeω(k)tq^0(kib)12iπD0+dkeω(k)tk(2ikb)ω(k)×[g0(0)+0tdseω(k)sg˙0(s)].

Proposition 2.2.

Assume that g˜0(t) is a smooth periodic function with period τ=2π and q(0,)g˜0S[0,). Suppose that g˜0(t) has the Fourier series given by

(27)g˜0(t)=n=aneint,t0,

where an with a0=0. Then, there is a unique periodic function g˜1(t) with period τ such that

(28)qx(0,t)g˜1(t)=O(t1),t,

and its Fourier series is given by

(29)g˜1(t)=n=(iknan)eint,

with

(30)kn={ib212κneiΘn/2,n>0,ib2+12κneiΘn/2,n<0,

where

κn=((b2+4a)2+16n2)14

and Θn=Arg(b24a4in) is the principal value of the argument.

Proof.

Let q(x,t) be the solution given in equation (25) with g0(t)=g˜0(t). We first prove equation (28). We write equation (26) as the sum of three integrals g1(t)=I1(t)+I2(t)+I3(t), where

I1(t)=12iπdkkeω(k)tq^0(k),I2(t)=12iπD+dkkeω(k)tq^0(kib),I3(t)=12iπD0+dkeω(k)tk(2ikb)ω(k)×[g0(0)+0tdseω(k)sg˙0(s)].

Note that using integration by parts, we find

I1(t)=12πdkeω(k)tq0(0)+12π0dxq˙0(x)dkeikxω(k)t=12πteat[eb2t4q0(0)+0dxe(x+bt)24tq˙0(x)],

which implies that I1(t)=O(eat) as t → . For the integral I2(t), we write I2(t) as

I2(t)=12iπ0dxq0(x)D0+dkke(ikb)xω(k)t.

By the Cauchy theorem, we deform the contour D0+ to (−, ) in the above equation, that is,

D0+dkke(ikb)xω(k)t=dkke(ikb)xω(k)t=πtexp[at(x+bt)24t].

Thus, we know that I2(t)=O(eat) as t → .

Regarding the integral I3(t), by integration by parts, we write I3(t) as

(31)I3(t)=12iπD0+dkk(2ikb)ω(k)×[g0(t)ω(k)0tdse(ts)ω(k)g0(s)].

Then, substituting the Fourier expansion (27) into (31), we find

(32)I3(t)=12iπD0+dkk(2ikb)×n0an(eintω(k)0tdse(ts)ω(k)+ins)=12iπD0+dkk(2ikb)×n0an(eintω(k)einteω(k)tω(k)+in).

For each n, kn defined in equation (30) solves the equation ω(k)+in=0 in D+. Thus, before we split the integral equation (32), we deform the contour D0+ to pass to the right of the singularities k = kn, denoted by D^0+. As a consequence, equation (32) can be written as

(33)I3(t)=12iπn0aneintD^0+dkk(2ikb)×(1ω(k)1ω(k)+in)12iπn0anD^0+dkeω(k)tk(2ikb)ω(k)+in.

Note that as k → 

k(2ikb)(1ω(k)1ω(k)+in)=2nk2+O(1k3).

Thus, by the residue theorem, we evaluate the first integral in equation (33) as

12iπn0aneintD^0+dkk(2ikb)×(1ω(k)1ω(k)+in)=n0iknaneint=g˜1(t).

For the second integral in equation (33), we use the steepest descent method for asymptotics of integrals. Note that −ω(k) has a critical point at k = −ib/2 and

k(2ikb)ω(k)+in=O(k+ib2),kib2.

Hence, we introduce the steepest descent contour CS, where CS is the line passing through the point k = −ib/2 and parallel to the real k-axis (depicted in Figure 1). Deforming the contour D^0+ to the contour CS and using the steepest descent method, we know that

D^0+dkeω(k)tk(2ikb)ω(k)+in=O(t1),t.

Therefore, g1(t)=g˜1(t)+O(t1) as t → . Moreover, the uniqueness of g˜1(t) follows from the fact that u(0,)g˜0(t)S[0,) (see [31] for some details). □

Remark 2.2.

It should be noted that the Fokas method can be applied effectively well to solve boundary problems with more general boundary conditions, such as the Neumann and Robin boundary conditions, namely, (i) Neumann boundary: qx(0,t)=g1(t); (ii) Robin boundary: αq(0,t)+qx(0,t)=h(t), α.

Neumann boundary condition. Solving equation (12) for g^0(k,t), the effective portion of the function g^0(k,t) is given by

(34)g^0(k,t)=1ik[q^0(kib)g^1(k,t)].

Thus, substituting equation (34) into equation (9), we find Q^(k,t) for equation (11) as

Q^(k,t)=1ik[(ikb)q^0(kib)(2ikb)g^1(k,t)].

Robin boundary condition. The Robin boundary value can be expressed in terms of the transforms

(35)αg^0(k,t)+g^1(k,t)=h^(k,t),

where h^(k,t)=0tdseω(k)sh(s). Substituting eqauation (12) into equation (35), the effective portion of g^0(k,t) is given by

g^0(k,t)=1ik+α[q^0(kib)h^(k,t)],

and then

g^1(k,t)=1ik+α(αq^0(kib)+ikh^(k,t)).

Using the above equations in equation (9), the effective portion for Q^(k,t) is given by

Q^(k,t)=1ik+α[(ikbα)q^0(kib)(2ikb)h^(k,t)].

3 The LAD equation on a finite interval

In Section 3, we discuss the Fokas method to solve the LAD equation posed on a finite interval, namely,

(36)qt=qxxbqxaq,0<x<L,t>0.

We denote the initial and boundary values as

q(x,0)=q0(x),q(0,t)=g0(t),q(L,t)=f0(t),qx(0,t)=g1(t),qx(L,t)=f1(t).

We assume that the boundary values g0(t) and f0(t) are given, while the boundary values g1(t) and f1(t) are unknown. We introduce the Fourier transform pair and the t-transform as

q^(k,t)=0Ldxeikxq(x,t),q(x,t)=12πdkeikxq^(k,t)

and

f^j(k,t)=0tdseω(k)sxjq(L,s),j=0,1.

Taking 0Ldx in equation (5) and employing 0tds in the resulting equation, equation (5) yields the following global relation

(37)eω(k)tq^(k,t)=q^0(k)G^(k,t)+eikLF^(k,t),

where

(38a)F^(k,t)=(ikb)f^0(k,t)+f^1(k,t),
(38b)G^(k,t)=(ikb)g^0(k,t)+g^1(k,t).

Applying the inverse Fourier transform, we find

(39)q(x,t)=12πdkeikxω(k)tq^0(k)12πdkeikxω(k)tG^(k,t)+12πdkeik(xL)ω(k)tF^(k,t).

By the Cauchy theorem, we deform the contour (−, ∞) to D+ for the integral involving the term G^(k,t) and to D for the integral involving the term F^(k,t) (noting x − L < 0), respectively. Then the reconstruction formula for q(x, t) is given by

(40)q(x,t)=12πdkeikxω(k)tq^0(k)12πD+dkeikxω(k)tG^(k,t)12πDdkeik(xL)ω(k)tF^(k,t).

Note that it requires to determine the unknown functions g^1(t) and f^1(t) in the representation of the solution equation (40). As before, the global relation equation (37) plays a crucial role. For simplicity, we write equation (37) as

(41)g^1(k,t)eikLf^1(k,t)=N(k,t)eω(k)tq^(k,t),

where

N(k,t)=q^0(k)(ikb)g^0(k,t)+eikL(ikb)f^0(k,t).

Replacing k → −k − ib in equation (41), we find

(42)g^1(k,t)e(ikb)Lf^1(k,t)=N(kib,t)eω(k)tq^(kib,t).

We solve the linear system of equations (41) and (42) for g^1(k,t) and f^1(k,t). As discussed in Section 2, note that the terms involving eω(k)tq^(k,t) and eω(k)tq^(kib,t) in equations (41) and (42) do not contribute in equation (40) thanks to the Cauchy theorem. Thus, the effective portions for the unknown functions g^1(k,t) and f^1(k,t) can be expressed as

(43)g^1(k,t)=1Δ(k)[eikLN(kib,t)e(ikb)LN(k,t)],
(44)f^1(k,t)=1Δ(k)[N(kib,t)N(k,t)],

where Δ(k)=eikLe(ikb)L is the determinant of the corresponding linear system. Finally, substituting equations (43) and (44) into equation (38), we find the functions G^(k,t) and F^(k,t) in equation (40) given by

(45)G^(k,t)=(ikb)g^0(k,t)+1Δ(k)×[eikLN(kib,t)e(ikb)LN(k,t)]F^(k,t)=(ikb)f^0(k,t)+1Δ(k)[N(kib,t)N(k,t)].

Remark 3.1.

Since the function Δ(k) has simple zeros at

km=mπLib2,m,

the functions F^(k,t) and G^(k,t) given in equation (45) have singularities at these points, which can be used to obtain an alternative series representation of the solution [16].

4 Examples

In Section 4, we consider the following Dirichlet boundary value problem

(46)qt=qxxqxq,x>0,t>0,

with the given initial and boundary values q(x,0)=q0(x) and q(0,t)=g0(t).

Example 4.1.

It is well-known that the function q(x,t)=etx solves equation (46). In this case,

q(x,0)=ex,q(0,t)=et.

We will derive the form of the solution directly from the above conditions. In terms of the transforms, we find

q^0(k)=1ik+1,g^0(k,t)=e(ω(k)+1)t1ω(k)+1,

where ω(k)=k2+ik+1. From equation (14), it follows that the solution is given by

(47)q(x,t)=12πD+dkeikxω(k)t(1ik+1+1ik2)12πD+dkeikx(2ik1)eteω(k)tω(k)+1,

where we have deformed the contour (−, ) to D+ by the Cauchy theorem. Note that the integrand of the first integral in equation (47) has a simple pole at k = i in D+, while the integrand of the second integral has a removable singularity at k = i in D+. Thus, by the residue theorem, we find q(x,t)=etx, as desired.

Example 4.2.

We consider the case

(48)q(x,0)=0,q(0,t)=C,

where C is a real constant. The function g^0(k,t) is given by

g^0(k,t)=C(eω(k)t1)ω(k),

which yields

(49)q(x,t)=C2πD+dkeikx(2ik1)1eω(k)tω(k).

The integral in equation (49) can not be computed explicitly, so we evaluate integral numerically. To this end, we use the efficient analytical-numerical scheme proposed in [34]. By the Cauchy theorem, we deform the contour D+ to L=L1L2 (depicted in Figure 1), where

L1={k=reiπ/6|r0},L2={k=re5iπ/6|r0}.

Note that along this contour L, the exponential terms in the integral decay rapidly for large k and hence numerical integration converges quickly. The accuracy and order of convergence for this hybrid analytical-numerical scheme has been discussed in [34]. Analytical and numerical solutions are shown in Figure 2 with C = 1. The analytical solution in Figure 2 (as blue solid curve) has been obtained by integrating equation (49) numerically along the contour L via the adaptive Gauss-Kronrod quadrature method with |k|=100 (or equivalently r = 100). The numerical solution displayed in Figure 2 (as red open dot) has been found by using MATLAB’s build-in function pdepe with x = 100. Figure 2 shows that the analytical solutions are in excellent agreement with numerical results.

Figure 2: Comparison of analytical solution (blue solid curve) and numerical solution (red open dot) with C = 1 and the initial and boundary values given in Equation (48). (a) Analytical and numerical solutions versus x when t = 0.2, 0.6 and 1.4 (from lower curve to upper one, respectively). (b) Analytical and numerical solutions versus t when x = 1, 2 and 4 (from upper curve to lower one, respectively) (color online).
Figure 2:

Comparison of analytical solution (blue solid curve) and numerical solution (red open dot) with C = 1 and the initial and boundary values given in Equation (48). (a) Analytical and numerical solutions versus x when t = 0.2, 0.6 and 1.4 (from lower curve to upper one, respectively). (b) Analytical and numerical solutions versus t when x = 1, 2 and 4 (from upper curve to lower one, respectively) (color online).

Example 4.3.

We consider the case of the periodic boundary value

(50)q(x,0)=0,q(0,t)=2sin2πt.

In this case, the solution is given by

(51)q(x,t)=12πD+dkeikxω(k)t(2ik1)g^0(k,t),

where

g^0(k,t)=2[2π+eω(k)t(ω(k)sin2πt2πcos2πt)]ω2(k)+4π2.

The integral in equation (51) can be computed by a similar way as discussed in Example 4.2. Figure 3a, b show that the analytical solution (shown in the figure as blue solid curve) is asymptotically periodic for large t and is again in excellent agreement with the numerical results (marked in the figure as red open dots).

Figure 3: Comparison of analytical solutions (blue solid curves) and numerical solutions (red open dots) versus t with the initial and boundary values given in equation (50), where x = 2 in (a) and x = 4 in (b), respectively (color online).
Figure 3:

Comparison of analytical solutions (blue solid curves) and numerical solutions (red open dots) versus t with the initial and boundary values given in equation (50), where x = 2 in (a) and x = 4 in (b), respectively (color online).

Example 4.4.

We consider the case that the Dirichlet boundary value is asymptotically periodic, namely,

(52)q(x,0)=0,q(0,t)=cos2πt1t4+1.

From equation (29), it follows that

(53)g˜1(t)=ik12e2iπt+ik12e2iπt,

where

k±1=i212(25+64π2)1/4ei2Θ±1

with Θ±1=Arg(58iπ). The asymptotic Neumann boundary value given in equation (53) is shown in Figure 4 by a blue solid curve. We have also compared the Neumann boundary value in this figure as a red dashed curve, which has been numerically obtained by the forward-difference derivative.

Figure 4: Comparison of numerically obtained Neumann boundary value (red dashed curve) and asymptotic formula (blue solid curve) given in equation (53) (color online).
Figure 4:

Comparison of numerically obtained Neumann boundary value (red dashed curve) and asymptotic formula (blue solid curve) given in equation (53) (color online).

Example 4.5.

We consider the LAD equation posed on the finite interval,

(54)qt=qxxqxq,0<x<1,t>0

with

q(x,0)=0,q(0,t)=2sin2πt,q(1,t)=0.

In this case, N(k,t)=(ik1)g^0(k,t), and hence

F^(k,t)=1Δ(k)(2ik1)g^0(k,t),G^(k,t)=1Δ(k)(2ik1)eikg^0(k,t),

where Δ(k)=eikeik1 and g^0(k,t) is given in Example 4.3. From equation (40), we find

q(x,t)=12πD+dkeik(x1)ω(k)t2ik1Δ(k)g^0(k,t)12πDdkeik(x1)ω(k)t2ik1Δ(k)g^0(k,t).

Using change of variable k → −k − i for the second integral, the solution q(x, t) can be written as

(55)q(x,t)=12πD+dkeω(k)t(eik(x1)e(ik1)(x1))×2ik1Δ(k)g^0(k,t).

The solution q(x, t) given in equation (55) is shown in Figure 5(a). We also have displayed the solution curves in Figure 5(b) as blue solid curves for x = 0.1, 0.5 and 0.8, which are in excellent agreement with the numerical solutions (red dashed curve).

Figure 5: (a) Analytical solution of equation (54) given in equation (55). (b) Comparison of the analytical solution (blue solid curve) and the numerical solution (red dashed curve) for x = 0.1, 0.5 and 0.8. (color online).
Figure 5:

(a) Analytical solution of equation (54) given in equation (55). (b) Comparison of the analytical solution (blue solid curve) and the numerical solution (red dashed curve) for x = 0.1, 0.5 and 0.8. (color online).

5 Concluding remarks

We have demonstrated the Fokas method to solve initial-boundary value problems for the LAD equation posed on the half line and a finite interval with general boundary conditions. In addition to solving the LAD equation, we have discussed the case of periodic boundary conditions, which commonly appears in physics and engineering. In particular, we have characterized the long time asymptotic behaviors of the solution and the unknown boundary value, which can be uniquely determined by the known asymptotically periodic boundary datum. These analytical predictions have been compared with numerical results showing the excellent agreement.

It should be noted that the presented method is relatively simple, but remarkably effective, to implement, finding an explicit integral representation of the solution. The method works equally well for boundary value problems for linear and nonlinear integrable systems. In contrast to classical transform methods such as the Fourier and Laplace transforms, the Fokas method has several advantages. The integral representation of the solution involves explicit (x, t)-dependence of the exponential form, and hence it allows to study the long time asymptotics of the solution. Also, from the integral representation of analytic functions, it makes possible to compute effectively the numerical solution [12], [20], [34]. More importantly, the method can be nonlinearizable in the sense that it is successfully used for analyzing nonlinear integrable PDEs as well as nonlinear lattices [35]. Recently, the method has been extended to study well-posedness, regularity and controllability of PDEs [36], [37], [38].

The LAD equation can be considered as a generalization of the heat equation and hence we expect that it could be possible to study boundary value problems with more complicated boundary data as shown in [39], [40], [41], [42], [43]. Also, note that the explicit representations of the solutions in finite and infinite domains can be used to analyze the effect of the exit and inlet boundary conditions for small Péclet number as discussed in [44]. We will address regarding these issues in the near future.


Corresponding author: Guenbo Hwang, Department of Mathematics and Institute of Natural Sciences, Daegu University, Gyeongsan, Gyeongbuk, 38453, Korea, E-mail:

Funding source: Daegu University

Award Identifier / Grant number: 2017

Acknowledgments

The author thanks anonymous referees for their useful comments and constructive suggestions. The work is supported by the Daegu University Research Grant 2017.

  1. Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.

  2. Research funding: This study was funded by Daegu University Research Grant 2017.

  3. Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

Appendix A

The inversion formula for the t-transform. The inversion formula for the spectral function equation (7) can be derived by performing spectral analysis of equation (4b). We consider the following spectral problem

(A.1)μt+ω(k)μ=g(t),0<t<T,

where μ=μ(t,k) and ω(k)=k2+ibk+a. Multiplying by the integrating factor eω(k)t, we found the Jost solutions given by

μ1(t,k)=0tdseω(k)(ts)g(s),μ2(t,k)=tTdseω(k)(ts)g(s).

Note that μ1 and μ2 are analytic and bounded in \D and D, respectively. Also, the jump condition, also known as a scalar Riemann-Hilbert problem, is given by

(A.2)(μ2μ1)(t,k)=eω(k)tg^(k,T),

where g^(k,t)=0tdseω(k)sg(s). By integration by parts, we find μ1,2=O(1/k2) as k → , and hence the solution of the Riemann-Hilbert problem equation (a.2) can be solved by the Plemelj formula [12]

μ(t,k)=12iπDdζeω(ζ)tg^(ζ,T)ζk,

where Ddζ=D+dζ+Ddζ. Substituting the above equation into equation (a.1), we find

g(t)=12iπDdζeω(ζ)tω(k)ω(ζ)ζkg^(ζ,T)=12iπDdζeω(ζ)t(ζ+ib+k)g^(ζ,T).

Letting k → −k − ib in the integral over D (and noting that if kD, then kibD+), we know that

Ddζeω(ζ)t(ζ+ib+k)g^(ζ,T)=D+dζeω(ζ)t(ζk)g^(ζ,T).

Thus, the reconstruction formula for g(t) can be found as

g(t)=12iπD+dζeω(ζ)t(2ζ+ib)g^(ζ,T).

Note that we can replace g^(ζ,T) by g^(ζ,t) in the above equation. Indeed, the difference between these terms

D+dζtTdseω(ζ)(ts)(2ζ+ib)g(s)

vanishes by the Cauchy theorem, since t − s < 0 and the integrand is analytic and bounded in D+. Therefore, we find the following inversion formula for the spectral function g^(k,t) as

(A.3)g(t)=12iπD+dkeω(k)tω(k)g^(k,t).

References

[1] M. T. van Genuchten and W. J. Alves, Analytical solutions of the one-dimensional convective-dispersive solute transport equation, Washington US Department of Agriculture, Technical Bulletins 157268, United States Department of Agriculture, Economic Research Service, 1982.Search in Google Scholar

[2] J. Zhong, C. Zeng, Y. Yuan, Y. Zhang, and Y. Zhang, “Numerical solution of the unsteady diffusion-convection-reaction equation based on improved spectral Galerkin method,” AIP Adv., vol. 8, p. 045314, 2018, https://doi.org/10.1063/1.5023332.Search in Google Scholar

[3] J. S. Pérez Guerrero, E. M. Pontedeiro, M. T. van Genuchten, and T. H. Skaggs, “Analytical solution of the one-dimensional advection–dispersion solute transport equation subject to time-dependent boundary conditions,” Chem. Eng. J., vol. 221, pp. 487–491, 2013, https://doi.org/10.1016/j.cej.2013.01.095.Search in Google Scholar

[4] I. Boztosun and A. Charafi, “An analysis of the linear advection–diffusion equation using mesh-free and mesh-dependent methods,” Eng. Anal. Bound. Elem., vol. 26, pp. 889–895, 2002, https://doi.org/10.1016/s0955-7997(02)00053-x.Search in Google Scholar

[5] A. Mojtabi and M. O. Deville, “One-dimensional linear advection–diffusion equation: analytical and finite element solutions,” Comput. Fluids, vol. 107, pp. 189–195, 2015, https://doi.org/10.1016/j.compfluid.2014.11.006.Search in Google Scholar

[6] A. Moranda, R. Cianci, and O. Paladino, “Analytical solutions of one-dimensional contaminant transport in soils with source production-decay,” Soil Syst., vol. 2, no. 40, pp. 1–16, 2018, https://doi.org/10.3390/soilsystems2030040.Search in Google Scholar

[7] D. H. Tang, E. O. Frind, and E. A. Sudicky, “Contaminant transport in fractured porous media: analytical solution for a single fracture,” Water Resour. Res., vol. 17, pp. 555–564, 1981, https://doi.org/10.1029/wr017i003p00555.Search in Google Scholar

[8] G. Dagan, “Theory of solute transport by groundwater,” Annu. Rev. Fluid Mech., vol. 19, pp. 183–215, 1987, https://doi.org/10.1146/annurev.fl.19.010187.001151.Search in Google Scholar

[9] J. S. Pérez Guerrero, L. C. G. Pimentel, T. H. Skaggs, and M. T. van Genuchten, “Analytical solution of the advection–diffusion transport equation using a change-of-variable and integral transform technique,” Int. J. Heat Mass Tran., vol. 52, pp. 3297–3304, 2009, https://doi.org/10.1016/j.ijheatmasstransfer.2009.02.002.Search in Google Scholar

[10] A. S. Fokas, “A unified transform method for solving linear and certain nonlinear PDEs,” Proc. Math. Phys. Eng. Sci., vol. 453, pp. 1411–1443, 1997, https://doi.org/10.1098/rspa.1997.0077.Search in Google Scholar

[11] A. S. Fokas, “Integrable nonlinear evolution equations on the half-line,” Commun. Math. Phys., vol. 230, pp. 1–39, 2002, https://doi.org/10.1007/s00220-002-0681-8.Search in Google Scholar

[12] A. S. Fokas, “A unified approach to boundary value problems,” in CBMS-NSF Regional Conf. Series in Applied Mathematics, Philadelphia, SIAM, 2008, https://doi.org/10.10.1137/1.9780898717068.10.1137/1.9780898717068Search in Google Scholar

[13] B. Deconinck, T. Trogdon, and V. Vasan, “The method of Fokas for solving linear partial differential equations,” SIAM Rev., vol. 56, pp. 159–186, 2014, https://doi.org/10.1137/110821871.Search in Google Scholar

[14] A. S. Fokas, “On the integrability of certain linear and nonlinear partial differential equations,” J. Math. Phys., vol. 41, pp. 4188–4237, 2000, https://doi.org/10.1063/1.533339.Search in Google Scholar

[15] A. S. Fokas, “A new transform method for evolution partial differential equations,” IMA J. Appl. Math., vol. 67, pp. 559–590, 2002, https://doi.org/10.1093/imamat/67.6.559.Search in Google Scholar

[16] A. S. Fokas and B. Pelloni, “A transform method for linear evolution PDEs on a finite interval,” IMA J. Appl. Math., vol. 70, pp. 564–587, 2005, https://doi.org/10.1093/imamat/hxh047.Search in Google Scholar

[17] M. J. Colbrook, N. Flyer, and B. Fornberg, “On the Fokas method for the solution of elliptic problems in both convex and non-convex polygonal domains,” J. Comput. Phys., vol. 374, pp. 996–1016, 2018, https://doi.org/10.1016/j.jcp.2018.08.005.Search in Google Scholar

[18] A. S. Fokas, “Two dimensional linear PDEs in a convex polygon,” Proc. Math. Phys. Eng. Sci., vol. 457, pp. 371–393, 2001, https://doi.org/10.1098/rspa.2000.0671.Search in Google Scholar

[19] B. Pelloni and D. A. Pinotsis, “The elliptic sine-Gordon equation in a half plane,” Nonlinearity, vol. 23, pp. 77–88, 2010, https://doi.org/10.1088/0951-7715/23/1/004.Search in Google Scholar

[20] M. J. Colbrook, A. S. Fokas, and P. Hashemzadeh, “A hybrid analytical-numerical technique for elliptic PDEs,” SIAM J. Sci. Comput., vol. 41, pp. A1066–A1090, 2019, https://doi.org/10.1137/18m1217309.Search in Google Scholar

[21] M. J. Colbrook, “Extending the unified transform: curvilinear polygons and variable coefficient PDEs,” IMA J. Numer. Anal., vol. 40, pp. 976–1004, 2020, https://doi.org/10.1093/imanum/dry085.Search in Google Scholar

[22] G. Biondini and G. Hwang, “Initial-boundary value problems for discrete evolution equations: discrete linear Schrödinger and integrable discrete nonlinear Schrödinger equations,” Inverse Probl., vol. 24, pp. 1–44, 2008,065011, https://doi.org/10.1088/0266-5611/24/6/065011.Search in Google Scholar

[23] B. Moon and G. Hwang, “Discrete linear evolution equations in a finite lattice,” J. Differ. Equ. Appl., vol. 25, pp. 630–646, 2019, https://doi.org/10.1080/10236198.2019.1613386.Search in Google Scholar

[24] A. S. Fokas, “Lax pair: a novel type of separability,” Inverse Probl., vol. 25, p. 123007, 2009, https://doi.org/10.1088/0266-5611/25/12/123007.Search in Google Scholar

[25] A. S. Fokas, “The generalized Dirichlet-to-Neumann map for certain nonlinear evolution PDEs,” Commun. Pure Appl. Math., vol. LVIII, pp. 639–670, 2005, https://doi.org/10.1002/cpa.20076.Search in Google Scholar

[26] D. G. Crowdy and E. Luca, “Solving Wiener-Hopf problems without kernel factorization,” Proc. Math. Phys. Eng. Sci., vol. 470, p. 20140304, 2014, https://doi.org/10.1098/rspa.2014.0304.Search in Google Scholar

[27] M. J. Colbrook, L. J. Ayton, and A. S. Fokas, “The unified transform for mixed boundary condition problems in unbounded domains,” Proc. R. Soc. A., vol. 475, p. 20180605, 2019, https://doi.org/10.1098/rspa.2018.0605.Search in Google Scholar

[28] A. B. de Monvel, A. Kotlyarov, D. Shepelsky, and C. Zheng, “Initial boundary value problems for integrable systems: towards the long-time asymptotics,” Nonlinearity, vol. 23, p. 2483, 2010, https://doi.org/10.1088/0951-7715/23/10/007.Search in Google Scholar

[29] J. Lenells and A. S. Fokas, “The unified method on the half-line: II. NLS on the half-line with t-periodic boundary conditions,” J. Phys. Math. Theor., vol. 45, p. 195202, 2012, https://doi.org/10.1088/1751-8113/45/19/195202.Search in Google Scholar

[30] G. Hwang and A. S. Fokas, “The modified Korteweg–de Vries equation on the half-line with a sine-wave as Dirichlet datum,” J. Nonlinear Math. Phys., vol. 20, pp. 135–157, 2013, https://doi.org/10.1080/14029251.2013.792492.Search in Google Scholar

[31] J. Lenells and A. S. Fokas, “The nonlinear Schrödinger equation with t-periodic data: II. Perturbative results,” Proc. Math. Phys. Eng. Sci., vol. 471, p. 20140926, 2015, https://doi.org/10.1098/rspa.2014.0926.Search in Google Scholar

[32] G. Hwang, “The modified Korteweg–de Vries equation on the quarter plane with t-periodic data,” J. Nonlinear Math. Phys., vol. 24, pp. 620–634, 2017, https://doi.org/10.1080/14029251.2017.1375695.Search in Google Scholar

[33] B. Moon and G. Hwang, “The Korteweg–de Vries equation on the quarter plane with asymptotically t-periodic data via the Fokas method,” Asymptot. Anal., vol. 107, pp. 115–133, 2018, https://doi.org/10.3233/asy-171452.Search in Google Scholar

[34] F. R. J. de Barros, M. J. Colbrook, and A. S. Fokas, “A hybrid analytical-numerical method for solving advection–dispersion problems on a half -line,” Int. J. Heat Mass Tran., vol. 139, pp. 482–491, 2019, https://doi.org/10.1016/j.ijheatmasstransfer.2019.05.018.Search in Google Scholar

[35] M. Toda, “Theory of a nonlinear lattice,” in Springer Series in Solid-State Science, vol. 20, Berlin, Springer-Verlag, 1989.10.1007/978-3-642-83219-2Search in Google Scholar

[36] A. S. Fokas, A. A. Himonas, and D. Mantzavinos, “The nonlinear Schrödinger equation on the half-line,” Trans. Amer. Math. Soc., vol. 369, pp. 681–709, 2017, https://doi.org/10.1090/tran/6734.Search in Google Scholar

[37] A. A. Himonas, D. Mantzavinos, and F. Yan, “The Korteweg–de Vries equation on an interval,” J. Math. Phys., vol. 60, p. 051507, 2019, https://doi.org/10.1063/1.5080366.Search in Google Scholar

[38] K. Kalimeris and T. Özsari, “An elementary proof of the lack of null controllability for the heat equation on the half line,” Appl. Math. Lett., vol. 104, p. 106241, 2020, https://doi.org/10.1016/j.aml.2020.106241.Search in Google Scholar

[39] D. Mantzavinos and A. S. Fokas, “The unified method for the heat equation: I. Non-separable boundary conditions and non-local constraints in one dimension,” Eur. J. Appl. Math., vol. 24, pp. 857–886, 2013, https://doi.org/10.1017/s0956792513000223.Search in Google Scholar

[40] A. S. Fokas and B Pelloni, “Generalized Dirichlet to Neumann map for moving initial-boundary value problems,” J. Math. Phys., vol. 48, p. 013502, 2007, https://doi.org/10.1063/1.2405405.Search in Google Scholar

[41] B. Pelloni and D. A. Smith, “Nonlocal and multipoint boundary value problems for linear evolution equations,” Stud. Appl. Math., vol. 141, pp. 46–88, 2018, https://doi.org/10.1111/sapm.12212.Search in Google Scholar

[42] B. Deconinck, B. Pelloni, and N. E. Sheils, “Non-steady-state heat conduction in composite walls,” Proc. Math. Phys. Eng. Sci., vol. 470, p. 20130605, 2014, https://doi.org/10.1098/rspa.2013.0605.Search in Google Scholar

[43] N. E. Sheils and B. Deconinck, “Heat conduction on the ring: interface problems with periodic boundary conditions,” App. Math. Lett., vol. 37, pp. 107–111, 2014, https://doi.org/10.1016/j.aml.2014.06.006.Search in Google Scholar

[44] M. Massabó, R. Cianci, and O. Paladino, “An analytical solution of the advection dispersion equation in a bounded domain and its application to laboratory experiments,” J. Appl. Math., vol. 2011, p. 493014, 2011, https://doi.org/10.1155/2011/493014.Search in Google Scholar

Received: 2020-04-13
Accepted: 2020-06-18
Published Online: 2020-07-31
Published in Print: 2020-08-27

© 2020 Walter de Gruyter GmbH, Berlin/Boston

Downloaded on 26.4.2024 from https://www.degruyter.com/document/doi/10.1515/zna-2020-0106/html
Scroll to top button