A Numerical Schemefor the Probability Density of the First Hitting Time for Some Random Processes

Abstract

Departing from a general stochastic model for a moving boundary problem, we consider the density function of probability for the first passing time. It is well known that the distribution of this random variable satisfies a problem ruled by an advection–diffusion system for which very few solutions are known in exact form. The model considers also a deterministic source, and the coefficients of this equation are functions with sufficient regularity. A numerical scheme is designed to estimate the solutions of the initial-boundary-value problem. We prove rigorously that the numerical model is capable of preserving the main characteristics of the solutions of the stochastic model, that is, positivity, boundedness and monotonicity. The scheme has spatial symmetry, and it is theoretically analyzed for consistency, stability and convergence. Some numerical simulations are carried out in this work to assess the capability of the discrete model to preserve the main structural features of the solutions of the model. Moreover, a numerical study confirms the efficiency of the scheme, in agreement with the mathematical results obtained in this work.

1. Preliminaries

In this manuscript, we let $B=\{B_s:s\ge 0\}$ be the Brownian motion, and we assume that $b,f,\sigma :[0,\infty )\to \mathbb{R}$ are functions. More concretely, we will suppose that f is continuous on $[0,\infty )$, $b\in L^1\left([0,\infty )\right)$ and $\sigma \in L^2\left([0,\infty )\right)$. The purpose of this work is to investigate the stochastic model

$$X_t^x=x+{\int }_0^{t}b\left(s\right)ds+{\int }_0^t\sigma \left(s\right)dB\left(s\right),\forall t\ge 0.$$

(1)

Here, x is some fixed (though arbitrary) real number. It is easy to check that the process $\{X_t^x:t\ge 0\}$ is Gaussian, and it has mean equal to ${\int }_0^{t}b\left(s\right)ds$ and variance ${\int }_0^t{\sigma }^2\left(s\right)ds$. Associated with this process, we define the first hitting time at the boundary f by

$${\tau }_x^f=inf\{t\ge 0:X_t^x=f\left(t\right)\},x<f\left(0\right).$$

(2)

It is indispensable to remember here that the variable ${\tau }_x^f$ has a continuous density if f has continuous derivatives up to the first order [1,2].

The class of stochastic processes considered in this investigation is a type of time-homogeneous Gauss-Markov diffusion process. As a consequence, the sample paths are continuous functions (which are nowhere differentiable) and, therefore, the concept of first hitting time through a boundary coincides with that of first passage time (see [3], for instance). In particular, this means that the random variable ${\tau }_x^f$ can be defined equivalently as

$${\tau }_x^f=inf\{t\ge 0:X_t^x\ge f\left(t\right)\},x<f\left(0\right).$$

(3)

Interestingly, the literature on the determination of first passage times in diffusion processes is vast. As examples, there are important research avenues which link the calculation of first passage time probability densities to the resolution of certain Volterra equations of the second kind for homogeneous [4] and non-homogeneous [5] diffusion processes that arise from Wiener processes.

The following result provides an analytical description of the probability density function of the first hitting time of the process $X_t^x$. In that result, we will employ the differential operator

$$A_t=\frac{{\left(x^2\sigma \left(t\right)\right)}^2}{2}\frac{{\partial }^2}{\partial x^2}+\left[(b\left(t\right)-f^{\prime }\left(t\right))x^2+x^3{\left(\sigma \left(t\right)\right)}^2\right]\frac{\partial }{\partial x}.$$

(4)

Theorem 1

(Macías-Díaz and Villa-Morales [6]). Let $u:\mathbb{R}\times [0,\infty )\to \mathbb{R}$ be bounded, and suppose that it satisfies the following parabolic problem with initial-boundary conditions:

$$\begin{array}{c}{\displaystyle \frac{\partial u}{\partial t}(x,t)=A_{t}u(x,t),\forall (x,t)\in \mathbb{R}\times {\mathbb{R}}^{+},}\\ \left\{\begin{array}{cc}u(x,0)=0,\hfill & \forall x\in \mathbb{R},\hfill \\ u(\infty ,t)=1,\hfill & \forall t>0.\hfill \end{array}\right.\end{array}$$

(5)

Then $P[{\tau }_x^f\le t]=u({(f\left(0\right)-x)}^{-1},t)$, for each $x<f\left(0\right)$ and $t>0$.

It is important to point out that, for the stochastic model considered in this work, the associated parabolic initial-boundary-value problem (5) is obtained from the Volterra integral equation of the second kind satisfied by the probability density of the first hitting time by making use of Itô’s formula and some properties on continuous martingales. The details on this relationship are provided in reference [6]. On the other hand, some exact solutions of this model are known in exact form [6]. More precisely, let $\alpha ,\beta \in \mathbb{R}$, $X_t=x+B_t$ and $f\left(t\right)=\alpha +\beta t$. For each $x<\alpha$,

$$P\left[{\tau }_x^f\le t\right]=1-\mathsf{\Phi }\left(\frac{\beta t+\alpha -x}{\sqrt{t}}\right)+e^{-2\left(\alpha -x\right)\beta }\mathsf{\Phi }\left(\frac{\beta t-\alpha +x}{\sqrt{t}}\right).$$

(6)

Here, $\mathsf{\Phi }$ is the cumulative distribution of a standard normal random variable. It is easy to check that

$$u(x,t)=1-\mathsf{\Phi }\left(\frac{\beta t+x^{-1}}{\sqrt{t}}\right)+e^{-2\beta x^{-1}}\mathsf{\Phi }\left(\frac{\beta t-x^{-1}}{\sqrt{t}}\right),\forall (x,t)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+},$$

(7)

solves (5) with $b\equiv 0$ and $\sigma \equiv 1$. However, for different expressions of the functions b, f and $\sigma$, the determination of explicit solutions for the problem (5) is a difficult task. This limitation motivates the design and analysis of numerical techniques to approximate the distribution of probability of the first hitting time in the stochastic model (1).

It is worth pointing out that the determination of the first hitting time often arises in problems associated with moving boundaries. Indeed, there are reports which study the first hitting time and the last exit time for Brownian motions to and from a moving boundary, respectively [7]. Other works tackle the investigation of the first hitting time for reducible diffusion processes [8], one-dimensional diffusion and compound Poisson processes [9], Markov processes between moving retaining or absorbing barriers [10], time-dependent Ornstein–Uhlenbeck processes to a moving boundary [11], moving boundaries for asymptotically stable Lévy processes [12], moving boundaries for some Gaussian stochastic processes [6], the pricing of continuously monitored barrier options under stochastic volatility with jumps [13], just to mention some examples [14,15,16].

In general, moving boundary problems have found interesting applications in the sciences and engineering. Indeed, this family of problems have been found in the diffusion of oxygen in absorbing tissues [17], in the analysis of progressive liquefaction [18], in the mathematical modeling of bread baking [19], in the asymptotic behavior of solutions of multidimensional problems associated with tumor growth [20], in the problem of thermal conduction in the solid phase of a phase-change material during melting driven by natural convection in the liquid [21] and problems on concrete carbonation [22]. Some recent studies have reported on applications of moving boundary problems to the radial fluid flow in infinite low-permeability reservoirs with threshold pressure gradients [23], to space-fractional diffusion logistic population models and density-dependent dispersal rates [24], to problems with variable specific heat and thermal conductivity [24], to the solvent diffusion within glassy polymers [25] and to American call options with transition costs [26], among other recent examples in the sciences and engineering [27,28].

In this manuscript, we investigate a moving boundary problem which considers the presence of a deterministic and a stochastic components. It is well known that the first hitting time has a distribution of probability governed by a parabolic advection–diffusion differential equation, for which very few solutions are known in exact form [6]. Motivated by this fact, we will propose a spatially symmetric finite-difference scheme to approximate the solutions of this differential model. Emphasis will be made on the preservation of the main analytical properties of the solutions of the problem, namely, the symmetry, the boundedness and the monotonicity. In that sense, the scheme proposed in this work will be a structure-preserving model [29,30]. Structure-preserving models have been designed to preserve the positivity and the symmetry of the solutions of Fisher-type equations [31], the monotonicity and boundedness of numerical model for Burgers–Huxley-type equations [32], the positivity of high-order Galerkin schemes for compressible Euler equations [33], the positivity and boundedness of numerical schemes for space-time fractional predator-prey models [34] and the energy and symmetry of Riesz space-fractional nonlinear wave equations [35]. In summary, the development of structure-preserving numerical techniques has been an important factor in problems where particular features of the solutions are physically meaningful [35,36].

In this work, we will investigate a generalization of the parabolic problem arising in Theorem 1. Precisely, we consider the generalized form of the parabolic equation of (5) given by the formula

$$\frac{\partial u(x,t)}{\partial t}=\frac{1}{2}\varphi (x,t)\frac{{\partial }^{2}u(x,t)}{\partial x^2}+\psi (x,t)\frac{\partial u(x,t)}{\partial x},\forall (x,t)\in \mathbb{R}\times {\mathbb{R}}^{+},$$

(8)

subjected to the initial and boundary conditions

$$\left\{\begin{array}{cc}u(x,0)=\chi \left(x\right),\hfill & \forall x\in \mathbb{R},\hfill \\ u(\infty ,t)=1,\hfill & \forall t>0.\hfill \end{array}\right.$$

(9)

We will assume that $\varphi ,\psi :\mathbb{R}\times [0,\infty )\to \mathbb{R}$ with $\varphi \ge 0$, and suppose that $\chi :\mathbb{R}\to \mathbb{R}$ is a nondecreasing function with the property that $0\le \chi \le 1$.

The present work is organized as follows. Section 2 is devoted to provide a structure-preserving numerical model to solve problem (13). Some equivalent presentations of the discrete system are provided, including a convenient vector form. We show therein that the discrete model conserves non-negativity, boundedness from above and monotonicity of the approximations. In Section 3 we establish that the numerical model is a consistent technique which is unconditionally stable and convergent. Some illustrative examples are provided in that section to show through computational simulations the most important features of the model. Finally, we close this manuscript with a section of discussions and some concluding remarks.

2. Discrete Model

For the remainder of this work, we will agree that $I_n=\{1,\dots ,n\}$ and ${\overline{I}}_n=I_n\cup \left\{0\right\}$, for each $n\in \mathbb{N}$. For computational purposes, we will let c, d and T be real numbers such that $c<d$ and $T>0$, and let $K,N\in \mathbb{N}$. We will fix a uniform partition of the interval $[c,d]$ consisting of K subintervals, of the form $c=x_0<x_1<\dots <x_k<\dots <x_K=d$, for each $k\in {\overline{I}}_K$. On the other hand, we will also fix a uniform partition of the temporal interval $[0,T]$ consisting of N subintervals, and we will let it take the form $0=t_0<t_1<\dots <t_n<\dots <t_N=T$, for each $n\in {\overline{I}}_N$. For each $(k,n)\in {\overline{I}}_K\times {\overline{I}}_N$, we let $u_k^n=u(x_k,t_n)$, and we use $U_k^n$ to represent a numerical approximation to the exact value of $u_k^n$.

Let $h=(d-c)/K$ and $\tau =T/N$. Following the usual conventions, we will employ ${\mathcal{R}}_h$ to represent the grid $\{x_k:k\in {\overline{I}}_K\}$. Moreover, ${\mathcal{V}}_h$ will denote the set of all real-valued functions defined on ${\mathcal{R}}_h$. It is easy to check then that $u^n={\left(u_k^n\right)}_{k=0}^K$ and $U^n={\left(U_k^n\right)}_{k=0}^K$ are members of ${\mathcal{V}}_h$, for each $n\in {\overline{I}}_N$. For the sake of convenience, we will let $u={\left(u^n\right)}_{n=0}^N$ and $U={\left(U^n\right)}_{n=0}^N$ in what follows. Moreover, we will require the following discrete linear operators.

Definition 1.

Let ${\left(V^n\right)}_{n=0}^N$ be a sequence in ${\mathcal{V}}_h$. We define the discrete linear operators

$$\begin{array}{cc}\hfill {\delta }_t{V}_k^n& =\frac{V_k^n-V_k^{n-1}}{\tau },\forall (k,n)\in {\overline{I}}_K\times I_{N-1},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\delta }_x^{\left(1\right)}{V}_k^n& =\frac{V_{k+1}-V_{k-1}}{2h},\forall (k,n)\in I_{K-1}\times {\overline{I}}_N,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\delta }_x^{\left(2\right)}{V}_k^n& =\frac{V_{k+1}-2V_k^n+V_{k-1}}{h^2},\forall (k,n)\in I_{K-1}\times {\overline{I}}_N.\hfill \end{array}$$

(12)

It is well known that these discrete operators provide consistent first-order or second-order approximations to suitable differential operators under appropriate regularity conditions. Moreover, with this nomenclature, the scheme to approximate the solutions of the initial-boundary-value problem (8) is given by the discrete system

$${\delta }_t{U}_k^{n+1}=\frac{1}{2}\varphi (x_k,t_{n+1}){\delta }_x^{\left(2\right)}{U}_k^{n+1}+\psi (x_k,t_{n+1}){\delta }_x^{\left(1\right)}{U}_k^{n+1},\forall (k,n)\in I_{K-1}\times {\overline{I}}_{N-1}.$$

(13)

About the boundary conditions, we will suppose that the solutions are prescribed at the boundary using Dirichlet data. More precisely, we will assume that ${\phi }_a,{\phi }_b:[0,T]\to \mathbb{R}$ are functions, and we will consider the following initial and boundary conditions:   

$$\left\{\begin{array}{cc}{U}_k^0=\chi \left(x_k\right),\hfill & \forall k\in I_{K-1},\hfill \\ U_0^n={\phi }_a\left(t_n\right),\hfill & \forall n\in {\overline{I}}_N,\hfill \\ U_K^n={\phi }_b\left(t_n\right),\hfill & \forall n\in {\overline{I}}_N.\hfill \end{array}\right.$$

(14)

For convenience, Figure 1 provides the forward-difference stencil of the discrete model (13) and (14). It is worth observing the spatial symmetry of the scheme at each time-step.

It is easy to check that the finite-difference scheme (13) and (14) is a two-step implicit technique. Moreover, the discrete model can be expressed in an equivalent form after some algebraic manipulations. For the sake of convenience, define the constants ${\varphi }_k^n=\varphi (x_k,t_n)$ and ${\psi }_k^n=\psi (x_k,t_n)$, for each $(k,n)\in {\overline{I}}_K\times {\overline{I}}_N$. Rearranging algebraically the terms in the expression of the general recursive formula in (13), we readily obtain the spatially symmetric implicit expression

$$-(R_k^{n+1}-r_k^{n+1})U_{k-1}^{n+1}+(1+2R_k^{n+1})U_k^{n+1}-(R_k^{n+1}+r_k^{n+1})U_{k+1}^{n+1}=U_k^n,\forall (k,n)\in I_{K-1}\times {\overline{I}}_{N-1}.$$

(15)

In this expression,

$$\begin{array}{cc}\hfill R_k^n& =\frac{\tau {\varphi }_k^n}{2h^2},\forall (k,n)\in I_{K-1}\times {\overline{I}}_N,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill r_k^n& =\frac{\tau {\psi }_k^n}{2h},\forall (k,n)\in I_{K-1}\times {\overline{I}}_N.\hfill \end{array}$$

(17)

It is obvious from (15) that the finite-difference scheme (13) is a linear model. In order to express it in vector form, we define de $(K+1)$-dimensional real vectors

$$\begin{array}{cc}\hfill {\mathbf{U}}^n& ={(U_0^n,U_1^n,\dots ,U_{K-1}^n,U_K^n)}^{\top },\forall n\in {\overline{I}}_N,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\mathbf{b}}^n& ={({\phi }_a\left(t_{n+1}\right),U_1^n,\dots ,U_{K-1}^n,{\phi }_b\left(t_{n+1}\right))}^{\top },\forall n\in {\overline{I}}_N.\hfill \end{array}$$

(19)

We have used ${}^{\top }$ here to represent the operation of matrix and vector transposition. It is obvious that ${\mathbf{U}}^n$ provides a vector representation of the function $U^n$, for each $n\in {\overline{I}}_N$. Moreover, notice that ${\mathbf{b}}^n$ includes already the data at the boundary, for each $n\in {\overline{I}}_N$. Now, for each $n\in {\overline{I}}_N$, we define the real tridiagonal matrix $A^n=\left(A_{ij}^n\right)$ of size $(K+1)\times (K+1)$ by

$$A^n=\left(\begin{array}{ccccccccc}1& 0& 0& 0& \cdots & 0& 0& 0& 0\\ -{\mu }_1^n& {\sigma }_1^n& -{\nu }_1^n& 0& \cdots & 0& 0& 0& 0\\ 0& -{\mu }_2^n& {\sigma }_2^n& -{\nu }_2^n& \cdots & 0& 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -{\mu }_{K-2}^n& -{\sigma }_{K-2}^n& -{\nu }_{K-2}^n& 0\\ 0& 0& 0& 0& \cdots & 0& -{\mu }_{K-1}^n& {\sigma }_{K-1}^n& -{\nu }_{K-1}^n\\ 0& 0& 0& 0& \cdots & 0& 0& 0& 1\end{array}\right).$$

(20)

Here, we let

$$\begin{array}{cc}\hfill {\mu }_k^n& =R_k^n-r_k^n,\forall (k,n)\in I_{K-1}\times {\overline{I}}_N,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\sigma }_k^n& =1+2R_k^n,\forall (k,n)\in I_{K-1}\times {\overline{I}}_N,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\nu }_k^n& =R_k^n+r_k^n,\forall (k,n)\in I_{K-1}\times {\overline{I}}_N.\hfill \end{array}$$

(23)

Using this nomenclature, it is easy to check that the spatially symmetric finite-difference scheme (13) and (14) can be equivalently rewritten in vector form as

$$\left\{\begin{array}{c}{A}^{n+1}{\mathbf{U}}^{n+1}={\mathbf{b}}^n,\forall n\in {\overline{I}}_{N-1},\hfill \\ \mathrm{such}\mathrm{that}{\mathbf{U}}^0=(\chi \left(x_0\right),\chi \left(x_1\right),\dots ,\chi \left(x_{K-1}\right),\chi \left(x_K\right)).\hfill \end{array}\right.$$

(24)

We will recall now some definitions from the literature.

Definition 2.

A square real matrix is an M-matrix if the following conditions are satisfied for the matrix:

(a)

it is strictly diagonally dominant,

(b)

its off-diagonal entries are less than or equal to zero, and 

(c)

its diagonal entries are greater than zero.

It is worth pointing out that every M-matrix is nonsingular, and that the components of its inverse are all positive numbers (see [37] and references therein). Now, we wish to establish conditions under which the matrices $A^n$ above are all M-matrices. To that end, we will show that the properties (a), (b) and (c) of Definition 2 are satisfied.

Lemma 1.

Let $n\in {\overline{I}}_N$. If $h{\psi }_k^n<{\varphi }_k^n$ is satisfied for each $k\in I_{K-1}$, then $A^n$ is an M-matrix.

Proof. 

Observe that the hypotheses imply that $R_k^n>r_k^n$ or, equivalently, that ${\mu }_k^n>0$, for each $k\in I_{K-1}$. Since ${\nu }_k^n>0$ also holds for each $k\in I_{K-1}$, it follows that the off-diagonal entries of $A^n$ are non-positive. It is obvious that the diagonal components of $A^n$ are positive, so it only remains to prove that the matrix is strictly diagonally dominant. Let $k\in I_{K-1}$, and consider the $(k+1)$th row of $A^n$. Note that

$$\sum _{\begin{array}{c}j=0\\ j\ne k\end{array}}^{K+1}|A_{kj}^n|=|{\mu }_k^n|+|{\nu }_k^n|=2R_k^n<1+2R_k^n={\sigma }_k^n=\left|A_{kk}^n\right|.$$

(25)

It readily follows then that $A^n$ is strictly diagonally dominant, so it is an M-matrix. □

Definition 3.

Let $\mathbf{V}$ and $\mathbf{W}$ be real vectors of the same dimension.

• We use the notation $\mathbf{V}\ge 0$ to represent that all the components of $\mathbf{V}$ are nonnegative.
• If the components of $\mathbf{V}$ are less than or equal to 1, then we denote this by $\mathbf{V}\le 1$. Note that $\mathbf{V}\le 1$ if and only if $\mathbf{1}-\mathbf{V}\ge 0$, where $\mathbf{1}$ is the vector of the same dimension of $\mathbf{V}$ whose components are equal to 1.
• If both $\mathbf{V}\ge 0$ and $\mathbf{V}\le 1$ are satisfied, then we will write $0\le \mathbf{V}\le 1$.
• We use the nomenclature $\mathbf{V}\le \mathbf{W}$ to denote that $\mathbf{W}-\mathbf{V}\ge 0$. Obviously, this represents that each component of $\mathbf{V}$ is less then or equal to the corresponding component of $\mathbf{W}$.
• Finally, we use $0\le \mathbf{V}\le \mathbf{W}\le 1$ to mean that $0\le \mathbf{V}$, $\mathbf{V}\le \mathbf{W}$ and $\mathbf{W}\le 1$ are all satisfied.

The following result establishes the existence and uniqueness of positive and bounded solutions of (13) and (14). The proof will hinge on various applications of Lemma 1.

Theorem 2.

Suppose that $0\le {\phi }_a\left(t\right)\le 1$ and $0\le {\phi }_b\left(t\right)\le 1$, for each $t\in [0,T]$, and let $0\le U^0\le 1$. Let $h{\psi }_k^n<{\varphi }_k^n$, for each $(k,n)\in I_{K-1}\times {\overline{I}}_N$. Then there exists a unique sequence ${\left(U^n\right)}_{n=0}^N$ in ${\mathcal{V}}_h$ satisfying (13), with the property that $0\le U^n\le 1$, for each $n\in {\overline{I}}_N$.

Proof. 

Beforehand, notice that the assumptions imply that the matrices $A^n$ are M-matrices, for each $n\in {\overline{I}}_N$. The proof of this result employs mathematical induction. Since the first approximation has the desired properties by assumption, we suppose that $0\le U^n\le 1$ holds for some $n\in {\overline{I}}_{N-1}$. The fact that $A^{n+1}$ is an M-matrix guarantees that it is nonsingular, and that the components of its inverse are all positive numbers. Moreover, the assumptions on the functions ${\phi }_a$ and ${\phi }_b$ and the induction hypothesis assure that ${\mathbf{b}}^n\ge 0$. As a consequence, ${\mathbf{U}}^{n+1}={\left(A^{n+1}\right)}^{-1}{\mathbf{b}}^n\ge 0$, and it only remains to prove that ${\mathbf{U}}^{n+1}\le 1$. To that end, let $\mathbf{1}$ be the $(K+1)$-dimensional vector all of whose components are equal to 1, and define ${\mathbf{V}}^{n+1}=\mathbf{1}-{\mathbf{U}}^{n+1}$ or, equivalently, ${\mathbf{U}}^{n+1}=\mathbf{1}-{\mathbf{V}}^{n+1}$. Substituting this identity into the vector equation of (24) and rearranging terms, we obtain that $A^{n+1}{\mathbf{V}}^{n+1}={\mathbf{c}}^{n+1}$, where

$${\mathbf{c}}^{n+1}=A^{n+1}\mathbf{1}-{\mathbf{b}}^n=\left(\begin{array}{c}1-{\phi }_a\left(t_{n+1}\right)\\ -{\mu }_1^{n+1}+{\sigma }_1^{n+1}-{\nu }_1^{n+1}-U_1^n\\ -{\mu }_2^{n+1}+{\sigma }_2^{n+1}-{\nu }_2^{n+1}-U_2^n\\ ⋮\\ -{\mu }_{K-1}^{n+1}+{\sigma }_{K-1}^{n+1}-{\nu }_{K-1}^{n+1}-U_{K-1}^n\\ 1-{\phi }_b\left(t_{n+1}\right)\end{array}\right)=\left(\begin{array}{c}1-{\phi }_a\left(t_{n+1}\right)\\ 1-U_1^n\\ 1-U_2^n\\ ⋮\\ 1-U_{K-1}^n\\ 1-{\phi }_b\left(t_{n+1}\right)\end{array}\right)\ge 0.$$

(26)

It readily follows that ${\mathbf{V}}^{n+1}\ge 0$ or, equivalently, that ${\mathbf{U}}^{n+1}\le 1$. Summarizing, we have checked that $0\le U^{n+1}\le 1$, and the result follows now by induction. □

The following result states that the finite-difference model (13) and (14) is a monotone technique.

Theorem 3.

Let $\chi ,\tilde{\chi }:[c,d]\to \mathbb{R}$ and ${\phi }_a,{\phi }_b,{\tilde{\phi }}_a,{\tilde{\phi }}_b:[0,T]\to \mathbb{R}$ be functions. Let $h{\psi }_k^n<{\varphi }_k^n$ and $h{\tilde{\psi }}_k^n<{\tilde{\varphi }}_k^n$, for each $(k,n)\in I_{K-1}\times {\overline{I}}_N$, and agree that ${\left(U^n\right)}_{n=0}^N$ and ${\left({\tilde{U}}^n\right)}_{n=0}^N$ are the solutions of the numerical model (13) and (14) corresponding to the initial-boundary data $(\chi ,{\phi }_a,{\phi }_b)$ and $(\tilde{\chi },{\tilde{\phi }}_a,{\tilde{\phi }}_b)$, respectively. Suppose that the initial-boundary conditions satisfy

(a)

$0\le \chi \left(x\right)\le \tilde{\chi }\left(x\right)\le 1$, for each $c\in [c,d]$, and 

(b)

$0\le {\phi }_a\left(t\right)\le {\tilde{\phi }}_a\left(t\right)\le 1$ and $0\le {\phi }_b\left(t\right)\le {\tilde{\phi }}_b\left(t\right)\le 1$, for each $t\in [0,T]$.

Then $0\le U^n\le {\tilde{U}}^n\le 1$, for each $n\in {\overline{I}}_N$.

Proof. 

Beforehand, observe that the hypotheses on the initial conditions assure that $0\le U^0\le 1$ and $0\le {\tilde{U}}^0\le 1$. Moreover, the boundary data are bounded in $[0,1]$. Theorem 2 implies now that the solutions ${\left(U^n\right)}_{n=0}^N$ and ${\left({\tilde{U}}^n\right)}_{n=0}^N$ of the problem (13) and (14) corresponding to the initial-boundary data $(\chi ,{\phi }_a,{\phi }_b)$ and $(\tilde{\chi },{\tilde{\phi }}_a,{\tilde{\phi }}_b)$ exist, are unique, and satisfy $0\le U^n\le 1$ and $0\le {\tilde{U}}^n\le 1$, for each $n\in {\overline{I}}_N$. It only remains to show that $U^n\le {\tilde{U}}^n$, for each $n\in {\overline{I}}_N$. Proceeding by induction, notice that this property is satisfied for $n=0$, so let us assume that it holds for some $n\in {\overline{I}}_{N-1}$. Observe that

$$\begin{array}{cc}\hfill A^{n+1}{\mathbf{U}}^{n+1}& ={\mathbf{b}}^n,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill A^{n+1}{\tilde{\mathbf{U}}}^{n+1}& ={\tilde{\mathbf{b}}}^n,\hfill \end{array}$$

(28)

are satisfied. Here,

$${\tilde{\mathbf{b}}}^n={({\tilde{\phi }}_a\left(t_{n+1}\right),{\tilde{U}}_1^n,\dots ,{\tilde{U}}_{K-1}^n,{\tilde{\phi }}_b\left(t_{n+1}\right))}^{\top }.$$

(29)

On the other hand, note that the induction hypothesis and the assumption on the boundary conditions yield that ${\tilde{\mathbf{b}}}^n-{\mathbf{b}}^n\ge 0$. Subtracting now (27) from (28), it follows that

$$A^{n+1}({\tilde{\mathbf{U}}}^{n+1}-{\mathbf{U}}^{n+1})={\tilde{\mathbf{b}}}^n-{\mathbf{b}}^n\ge 0.$$

(30)

We use the fact now that $A^{n+1}$ is an M-matrix under the hypotheses of this theorem. The fact that the components of the inverse of $A^{n+1}$ are all positive numbers establishes that ${\tilde{\mathbf{U}}}^{n+1}-{\mathbf{U}}^{n+1}\ge 0$ or, equivalently, that ${\tilde{\mathbf{U}}}^{n+1}\ge {\mathbf{U}}^{n+1}$, as desired. □

The following is a straightforward consequence of Theorem 3. Its relevance arises from the fact that the solutions of (5) are cumulative distributions of probability in the sense of Theorem 1.

Corollary 1.

Suppose that $h{\psi }_k^n<{\varphi }_k^n$ for each $(k,n)\in I_{K-1}\times {\overline{I}}_N$. Assume that $\chi :[c,d]\to \mathbb{R}$, and let ${\phi }_a,{\phi }_b:[0,T]\to \mathbb{R}$ be nondecreasing functions satisfying

(a)

$0\le \chi \left(x\right)\le 1$, for each $c\in [c,d]$, and 

(b)

$0\le {\phi }_a\left(t\right)\le 1$ and $0\le {\phi }_b\left(t\right)\le 1$, for each $t\in [0,T]$.

If $U^0\le U^1$ then the solution ${\left(U^n\right)}_{n=0}^N$ of (13) and (14) satisfies $0\le U^n\le U^{n+1}\le 1$, for each $n\in {\overline{I}}_{N-1}$.

Proof. 

To start with, notice that the existence and the uniqueness of nonnegative and bounded solutions is guaranteed by Theorem 2, so it only remains to show that $U^n\le U^{n+1}$, for each $n\in {\overline{I}}_{N-1}$. To that end, let ${\tilde{U}}^n=U^{n+1}$, for each $n\in {\overline{I}}_{N-1}$. An application of Theorem 3 to the sequences ${\left(U^n\right)}_{n=0}^{N-1}$ and ${\left({\tilde{U}}^n\right)}_{n=0}^{N-1}$ establishes now that $U^n\le {\tilde{U}}^n$ is satisfied for each $n\in {\overline{I}}_{N-1}$, whence the conclusion of this result readily follows. □

3. Numerical Properties

The present section is devoted to study the main numerical features of the finite-difference scheme (13) and (14). In particular, we will investigate the consistency of the scheme, its stability and its convergence. In a first stage, we will tackle the problem of the consistency of our numerical technique. To that end, let $\mathsf{\Omega }=(c,d)\times (0,T)$ and define the differential operator

$$\mathcal{L}\left(u_k^n\right)=\frac{\partial u(x_k,t_n)}{\partial t}-\frac{1}{2}\varphi (x_k,t_n)\frac{{\partial }^{2}u(x_k,t_n)}{\partial x^2}-\psi (x_k,t_n)\frac{\partial u(x_k,t_n)}{\partial x},\forall (k,n)\in I_{K-1}\times I_N,$$

(31)

and the difference operator

$$L\left(u_k^n\right)={\delta }_t{u}_k^n-\frac{1}{2}\varphi (x_k,t_n){\delta }_x^{\left(2\right)}{u}_k^n-\psi (x_k,t_n){\delta }_x^{\left(1\right)}{u}_k^n,\forall (k,n)\in I_{K-1}\times I_N.$$

(32)

For the sake of convenience, we let $\mathcal{L}\left(u^n\right)={\left(\mathcal{L}\left(u_k^n\right)\right)}_{k=0}^K$ and $L\left(u^n\right)={\left(L\left(u_k^n\right)\right)}_{k=0}^K$, for each $n\in I_N$.

Definition 4.

Define the function ${\parallel ·\parallel }_{\infty }:{\mathcal{V}}_h\to \mathbb{R}$ by ${\parallel V\parallel }_{\infty }=max\left\{\right|V_k|:k\in I_{K-1}\}$, for each $V\in {\mathcal{V}}_h$.

Theorem 4.

Let $u\in {\mathcal{C}}_{x,t}^{4,2}\left(\overline{\mathsf{\Omega }}\right)$, and suppose that ϕ and ψ are bounded in $\overline{\mathsf{\Omega }}$. There exists a constant $C\ge 0$ which is independent of τ and h, such that $\parallel \mathcal{L}\left(u^n\right)-L\left(u^n\right){\parallel }_{\infty }\le C(\tau +h^2)$, for each $n\in I_N$.

Proof. 

Using the regularity of the functions u, $\varphi$ and $\psi$ along with Taylor’s theorem, it is easy to show that there exist nonnegative constants $C_1$, $C_2$ and $C_3$ which are independent of $\tau$ and h, such that

$$\begin{array}{cc}\hfill \left|\frac{\partial u(x_k,t_n)}{\partial t}-{\delta }_t{u}_k^n\right|& \le C_1\tau ,\forall (k,n)\in I_{K-1}\times I_N,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \left|\varphi (x_k,t_n)\frac{{\partial }^{2}u}{\partial x^2}(x_k,t_n)-\varphi (x_k,t_n){\delta }_x^{\left(2\right)}{u}_k^n\right|& \le C_2{h}^2,\forall (k,n)\in I_{K-1}\times I_N,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \left|\psi (x_k,t_n)\frac{\partial u}{\partial x}(x_k,t_n)-\psi (x_k,t_n){\delta }_x^{\left(1\right)}{u}_k^n\right|& \le C_3{h}^2,\forall (k,n)\in I_{K-1}\times I_N.\hfill \end{array}$$

(35)

Using the triangle inequality, it is easy to check that

$$\left|\mathcal{L}\left(u_k^n\right)-L\left(u_k^n\right)\right|\le C_1\tau +\frac{1}{2}{C}_2{h}^2+C_3{h}^2\le C(\tau +h^2),\forall (k,n)\in I_{K-1}\times I_N.$$

(36)

The conclusion of this theorem readily follows now if we define $C=C_1+\frac{1}{2}{C}_2+C_3$, which is a nonnegative constant that is independent of $\tau$ and h. □

Next, we tackle the problems on the stability and the convergence of the finite-difference scheme. The following technical results will be of utmost importance to that end.

Lemma 2

(Chen et al. [38]). Suppose that A is a real matrix of size $(K+1)\times (K+1)$ which satisfies

$$\sum _{\begin{array}{c}j=1\\ j\ne k\end{array}}^{K+1}|a_{kj}|\le |a_{kk}|-1,\forall k\in I_{K+1}.$$

(37)

Then ${\parallel \mathbf{V}\parallel }_{\infty }\le {∥A\mathbf{V}∥}_{\infty }$ is satisfied for all $\mathbf{V}\in {\mathbb{R}}^{K+1}$.

Lemma 3.

Let $n\in {\overline{I}}_N$. If $h{\psi }_k^n<{\varphi }_k^n$, for each $k\in I_{K-1}$, then ${\parallel \mathbf{V}\parallel }_{\infty }\le {∥A^n\mathbf{V}∥}_{\infty }$, for each $\mathbf{V}\in {\mathbb{R}}^m$.

Proof. 

Let $k\in I_{K+1}$. Observe that when $2\le k\le K$,

$$\sum _{\begin{array}{c}j=1\\ j\ne k\end{array}}^{K+1}|A_{kj}^n|=2R_k^n=\left|A_{kk}^n\right|-1,$$

(38)

which means that the inequality (37) is satisfied in this case. It is easy to check that this inequality is also satisfied when $k=1$ or $k=K+1$. The conclusion follows now from Lemma 2.

We establish now the nonlinear stability of the finite-difference model (13) and (14). To that end, we will consider two sets of initial-boundary conditions, which will be denoted by $(\chi ,{\phi }_a,{\phi }_b)$ and $(\tilde{\chi },{\tilde{\phi }}_a,{\tilde{\phi }}_b)$. The corresponding solutions obtained using the numerical model will be denoted by ${\left({\mathbf{U}}^n\right)}_{n=0}^N$ and ${\left({\tilde{\mathbf{U}}}^n\right)}_{n=0}^N$, respectively.

Theorem 5.

Suppose that $h{\psi }_k^n<{\varphi }_k^n$ holds for each $(k,n)\in I_{K-1}\times {\overline{I}}_N$, and let ${\left(U^n\right)}_{n=0}^N$ and ${\left({\tilde{U}}^n\right)}_{n=0}^N$ be solutions of the discrete model (13) and (14) corresponding to the initial-boundary conditions $(\chi ,{\phi }_a,{\phi }_b)$ and $(\tilde{\chi },{\tilde{\phi }}_a,{\tilde{\phi }}_b)$, respectively. Then $\parallel {\tilde{U}}^n-U^n{\parallel }_{\infty }\le {\parallel {\tilde{U}}^0-U^0\parallel }_{\infty }$, for each $n\in {\overline{I}}_N$.

Proof. 

Under the conditions of this result, the conclusion of Lemma 3 is satisfied. As a consequence

$$\parallel {\tilde{\mathbf{U}}}^{n+1}-{\mathbf{U}}^{n+1}{\parallel }_{\infty }\le \parallel A^{n+1}({\tilde{\mathbf{U}}}^{n+1}-{\mathbf{U}}^{n+1}){\parallel }_{\infty }=\parallel {\tilde{\mathbf{b}}}^n-{\mathbf{b}}^n{\parallel }_{\infty }={\parallel {\tilde{\mathbf{U}}}^n-{\mathbf{U}}^n\parallel }_{\infty },\forall n\in {\overline{I}}_{N-1}.$$

(39)

The conclusion of this result readily follows now using mathematical induction. □

The following result will be required to prove the linear stability of our scheme. In that result, the function $\rho (·)$ will denote the spectral radius of real and square matrices.

Lemma 4

(Tian and Huang [39]). Let $M=\left(m_{ij}\right)$ be an M-matrix, and let $N=\left(n_{ij}\right)$ be a nonnegative matrix of the same size as M. If M is strictly diagonally dominant by rows then $\rho \left(M^{-1}N\right)$ satisfies

$$\begin{array}{c}\hfill {\displaystyle \rho \left(M^{-1}N\right)\le \underset{i\in N}{max}\left\{\frac{{\sum }_{j=1}^n{n}_{ij}}{m_{ii}+{\sum }_{j\ne i}{m}_{ij}}\right\}.}\end{array}$$

(40)

Theorem 6.

If $h{\psi }_k^n<{\varphi }_k^n$ is satisfied for each $(k,n)\in I_{K-1}\times {\overline{I}}_N$, then the scheme (13) and (14) is linearly stable.

Proof. 

The hypotheses guarantee that $A^n$ is an M-matrix, for each $n\in {\overline{I}}_N$. On the other hand, the applying the previous lemma with N being the identity matrix of the same size as $A^n$ yields

$$\rho \left({\left(A^n\right)}^{-1}\right)\le \underset{i\in I_{K+1}}{max}\left\{\frac{1}{1+2R_i^n-(R_i^n-r_i^n)-(R_i^n+r_i^n)}\right\}=1,\forall n\in {\overline{I}}_N.$$

(41)

The linear stability of the scheme readily follows from this fact. □

The next result is a direct consequence of the consistency and stability properties of the scheme.

Theorem 7.

Let $u\in {\mathcal{C}}_{x,t}^{4,2}\left(\overline{\mathsf{\Omega }}\right)$ be a solution of the problem (8) and (9), and suppose that ϕ and ψ are bounded in $\overline{\mathsf{\Omega }}$. The solutions of the corresponding discrete model (13) and (14) converge to u with order $\mathcal{O}(\tau +h^2)$, whenever the parameter h satisfies $h{\psi }_k^n<{\varphi }_k^n$, for each $(k,n)\in I_{K-1}\times {\overline{I}}_N$. □

Before closing this section, we will provide some illustrative simulations on the performance of the numerical model proposed in this work. In particular, we wish to exhibit the capability of the numerical model to preserve the positivity, the boundedness and the monotonicity of the discrete solutions. To that end, we will compare the numerical results against the exact solution described after Theorem 1. It is worth pointing out that the initial and the boundary data will be prescribed exactly using that analytical solution.

In our simulations, we will consider the stochastic process $X_t=x+B_t$ with the moving boundary $f\left(t\right)=\alpha +\beta t$. In that case, the solution of the problem (5) is given by the function (7). We will set now $\alpha =1$, $\mathsf{\Omega }=(0,1)\times (0,10)$, and employ the computational parameters $h=0.01$ and $\tau =0.0001$. Under these circumstances, Figure 2 shows the exact solution (left) and the numerical approximation (right) to the solution of this problem. We used the finite-difference scheme (13) and (14) together with the parameter value $\beta =0.25$ (top row), $\beta =0.5$ (middle row) and $\beta =1$ (bottom row). The results show a good qualitative agreement between the numerical and the analytical solutions of this problem over the space-time domain. For comparison purposes, we have plotted the exact and numerical solutions of this problem at time $T=10$ using (a) $\beta =0.25$, (b) $\beta =0.5$, (c) $\beta =0.75$ and (d) $\beta =1$. The solutions are plotted in Figure 3, and they prove that there exists a good agreement between the numerical and the theoretical results, as expected. To check this fact, Figure 4 provides absolute error plots corresponding to the graphs in Figure 3. Finally, Figure 5 provides a graph of the analysis of precision (using the ${\parallel ·\parallel }_{\infty }$ norm) versus computational cost, and contrasts the results against those obtained in [6] for the model described in the caption therein. The present scheme shows a superior performance due to the fact that less computer operations are required.

4. Discussion and Conclusions

To start with, we present various remarks with the intention of responding various important criticisms and suggestions by the reviewers of this work.

• Regarding the calculation of the first hitting time densities, only some few solutions are known in exact form for certain types of boundaries [4,5]. In every other case, numerical schemes should be considered and, in that sense, the present work contributes to the development of reliable numerical models to approximate the densities of first hitting times. However, we must point out that those densities can also be approximated by numerical schemes derived from the approximation to the associated integral equations. It is worth pointing out that some software has been developed following that approach [40]. The author of the present manuscript has followed a different approach, but the comparison between the present scheme and that reported in [40] is an interesting topic of investigation. However, extensive tests would need to be carried out under similar computational conditions. To start with, the codes must be in the same language and tests must be performed under the same computational architecture. To this day, the author has not devoted time to perform these tests, but it is an interesting task that the author of this work intends to carry out in the future.
• We must mention that it is possible to generalize the numerical procedure proposed in this work to account for more general stochastic processes than that considered herein, though some suitable adjustments are needed to that effect. As an example, it is well known that the probability distribution of the first passage time for Lévy time-charged Brownian motion with drift is a solution of a time-fractional advection–diffusion equation subjected to initial-boundary conditions [41]. In that case, the present approach can be modified to account for the presence of Caputo time-fractional derivatives and approximate the solutions of the associated partial differential equation. However, the structural and numerical analysis of that scheme would require a substantial amount of more effort and time. That could be the research topic of some independent report in the near future.
• It is important to point out that the same mathematical problem was considered in the work [6]. In that paper, the authors proposed a finite-difference scheme to approximate the probability distribution of the first hitting time associated with the same problem investigated in the present manuscript. The present approach has more advantages, though. To start with, the current numerical model is easier to implement computationally than that investigated in [6]. This fact can be easily checked from the Matlab code in Appendix A. Moreover, various theoretical and numerical advantages in the present approach are described below in the conclusions.
• The numerical scheme presented in this work may be relatively simple, but the analysis is rigorous and has various advantages. For example, it preserves the positivity, the boundedness and the monotonicity of the numerical solutions. This is an important feature in view that the relevant solutions of the parabolic problem considered herein possess these properties. The existence and uniqueness of solutions was thoroughly established here. Moreover, the scheme proposed in this work was theoretically analyzed for the numerical properties. Another advantage of our present approach is that the proofs of stability and convergence do not require additional conditions other than those asked to guarantee the existence of solutions.

Summarizing, we proposed a finite-difference model to solve a parabolic initial-value-problem associated with a stochastic model. The motivating problem is a stochastic equation which considers both the presence of a deterministic component and a Gaussian term, as well as nonconstant coefficients. This system has various applications, including the modeling of moving boundary systems and crossing times. It is known that the first hitting time of the stochastic model has a probability distribution which is described by a linear parabolic initial-boundary-value problem. The exact solution of that problem is only known in some specific scenarios, whence the need to propose a reliable numerical model to solve it is an important problem of investigation. Moreover, in view that the relevant solutions of the initial-boundary-value problem are cumulative distributions of probability, it is important to guarantee that the numerical solutions be nondecreasing and nonnegative functions, bounded from above by 1.

The proposed finite-difference scheme is a numerical model which is capable of preserving the properties of monotonicity, nonnegativity and boundedness, and that has spatial symmetry. These features are thoroughly established in this work using matrix theory. Solutions of the discrete model exist, and the computer implementation of this discrete model can be carried out efficiently using matrix representations. The scheme is quadratically consistent in space, and linearly consistent in time. The properties of linear stability, nonlinear stability and convergence of the scheme are also theoretically elucidated. In particular, we show that our numerical model has order of convergence equal to 2 in space, and equal to 1 in time. We implemented the scheme computationally, and some computer simulations were provided in this work in order to show that the scheme is capable of preserving the main features of the relevant solutions of the initial-boundary-value problem. In addition, our computational results confirm that the scheme is stable and convergent.

The present approach has various advantages with respect to some previous efforts reported in the literature to approximate the cumulative probability distribution of the first hitting time to similar moving boundary problems. For example, in reference [6], the authors proposed and analyzed a finite-difference scheme to approximate the probability distribution of the first hitting time for the same problem investigated in the present work. Various advantages are found in the present manuscript when compared against that work of the literature:

• The numerical model introduced in the present work has a much easier computational implementation, in view that the system to solve at each time step is simpler. Indeed, the simplified Matlab program to obtain the simulations in this manuscript is that presented in Appendix A. Notice the simplicity of the code. Obviously, the computer program may be simplified or generalized even more.
• The current discretization presents spatial symmetry.
• The analysis of stability is the numerical model of the present manuscript was carried out without assuming that the matrices $A^{n+1}$ depend on the approximations at the time $t_n$. This is a strong limitation of the results in [6]. Indeed, in that work, the authors needed to assume that the matrices did not depend on ${\mathbf{U}}^n$, which was an unrealistic assumption.
• The limitations of [6] described in the previous remark carried over to the analysis of convergence. Indeed, the model proposed in that paper had order of convergence $\mathcal{O}({\tau }^2+h^2)$, but the hypotheses to guarantee that order were utterly restrictive and unrealistic.
• Stability and convergence were established in the present discretization. On the contrary, the paper mentioned above requires additional constraints which limit the approach severely. Our present approach has convergence order $\mathcal{O}(\tau +h^2)$, but the theoretical result need less constraints to guarantee this order. We have performed numerical tests of convergence (not presented here in order to avoid redundancy), and they show that the present scheme has linear order of convergence in the temporal variable, and second order in the spatial domain.
• It is worth mentioning that there are various reports available in the literature in which problems with moving boundaries have been investigated numerically [42,43,44,45,46,47,48]. However, the present is an easy-to-implement model which preserves various advantages from the structural and numerical points of view.