A discontinuous Galerkin method for approximating the stationary distribution of stochastic fluid-fluid processes

TL;DR

This paper develops a discontinuous Galerkin numerical method to approximate the stationary distribution of stochastic fluid-fluid processes, enabling practical computation where previous formulas were analytically explicit but numerically intractable.

Contribution

It introduces a novel discontinuous Galerkin approach for numerical approximation of the stationary distribution in stochastic fluid-fluid processes.

Findings

Successfully applied to an on-off bandwidth sharing system

Demonstrates accurate approximation of stationary distributions

Provides a practical computational tool for complex stochastic processes

Abstract

Introduced by Bean and O'Reilly (2014), a stochastic fluid-fluid process is a Markov processes $\{X_t, Y_t, \varphi_t\}_{t \geq 0}$, where the first fluid $X_t$ is driven by the Markov chain $\varphi_t$, and the second fluid $Y_t$ is driven by $\varphi_t$ as well as by $X_t$. That paper derived a closed-form expression for the joint stationary distribution, given in terms of operators acting on measures, which does not lend itself easily to numerical computations. Here, we construct a discontinuous Galerkin method for approximating this stationary distribution, and illustrate the methodology using an on-off bandwidth sharing system, which is a special case of a stochastic fluid-fluid process.

Tables

Table 1. Table 1 : Output rates ζ 1 subscript 𝜁 1 \zeta_{1} and ζ 2 subscript 𝜁 2 \zeta_{2} for the buffers in different scenarios, as specified in the model in [ 6 ] . Note that while Buffer 1 is independent of Buffer 2, its output rate ζ 1 subscript 𝜁 1 \zeta_{1} depends on X t subscript 𝑋 𝑡 X_{t} .

Buffer $1$	Buffer $2$ and Phase	${\zeta }_1$	${\zeta }_2$
$X_1>x^{\ast }$	$\{Y_t\in [0,\infty )$, ${\phi }_t\in 𝒮\}$	$\kappa$	$0$
	$\left\{Y_t>0,{\phi }_t\in 𝒮\right\}$ or $\left\{Y_t=0,{\phi }_t\in \{11,01\}\right\}$	${\theta }_1$	${\theta }_2$
	$\left\{Y_t=0,{\phi }_t\in \{10,00\}\right\}$	${\theta }_1$	$0$
$X_t=0$	$\left\{Y_t>0,{\phi }_t\in 𝒮\right\}$ or $\left\{Y_t=0,{\phi }_t\in \{11,01\}\right\}$	0	$\kappa$
	$\left\{Y_t=0,{\phi }_t\in \{10,00\}\right\}$	$0$	$0$

Table 2. Table 2 : Output rates ζ 1 subscript 𝜁 1 \zeta_{1} and ζ 2 subscript 𝜁 2 \zeta_{2} for the buffers in different scenarios, as modified slightly from the model in [ 6 ] . Note that the output rate ζ 1 subscript 𝜁 1 \zeta_{1} no longer depends on the value of X t subscript 𝑋 𝑡 X_{t} , except for a boundary condition at 0 0 .

Buffer $1$	Buffer $2$ and Phase	${\zeta }_1$	${\zeta }_2$
$X_1>x^{\ast }$	$\{Y_t\in [0,\infty )$, ${\phi }_t\in 𝒮\}$	${\theta }_1$	$0$
	$\left\{Y_t>0,{\phi }_t\in 𝒮\right\}$ or $\left\{Y_t=0,{\phi }_t\in \{11,01\}\right\}$	${\theta }_1$	${\theta }_2$
	$\left\{Y_t=0,{\phi }_t\in \{10,00\}\right\}$	${\theta }_1$	$0$
$X_t=0$	$\left\{Y_t>0,{\phi }_t\in 𝒮\right\}$ or $\left\{Y_t=0,{\phi }_t\in \{11,01\}\right\}$	0	$\kappa$
	$\left\{Y_t=0,{\phi }_t\in \{10,00\}\right\}$	$0$	$0$

Table 3. Table 3 : The functions χ ^ 0 ​ ( x ) superscript ^ 𝜒 0 𝑥 \widehat{\chi}^{0}(x) and χ ^ + superscript ^ 𝜒 \widehat{\chi}^{+} are piecewise linear DG approximations of the limiting marginal densities χ 0 ​ ( x ) superscript 𝜒 0 𝑥 \chi^{0}(x) and χ + ​ ( x ) superscript 𝜒 𝑥 \chi^{+}(x) over the stencil 𝝎 43 , 0.4 , 0.001 . subscript 𝝎 43 0.4 0.001 \bm{\omega}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}43},0.4,0.001}.

		${\alpha }_2=11,{\beta }_2=1$	${\alpha }_2=16,{\beta }_2=1$	${\alpha }_2=22,{\beta }_2=1$
${\int }_{[0,ℐ]}{\hat{\chi }}^0\left(x\right)dx$		$\approx$ 0.0	0.184	0.312
${\int }_{[0,ℐ]}{\hat{\chi }}^{+}\left(x\right)dx$		$\approx$ 1.0	0.816	0.688

Table 4. Table 4 : Computational times and storage comparisons between piecewise-linear and -constant approximations. Overall storage is the total storage of all the operators from ( 7 ) to ( 24 ). The computations were performed on 2.5Ghz Intel Core i7 with 16GB of RAM running OSX 10.10.5. The code was implemented in python, using scientific python libraries.

Basis	$h$	Error	Comp. Times	# of Elements	Overall Storage
Piecewise Constant	1.5	0.58	0.21 sec	88	0.5 MB
	0.5	0.25	0.31 sec	248	4.1 MB
	0.05	0.03	23 sec	2408	380 MB
Piecewise Linear	1.5	0.01	0.21 sec	168	2.1 MB
	0.5	$1.1\times {10}^{-3}$	0.78 sec	488	16 MB
	0.05	$2.1\times {10}^{-5}$	130 sec	4816	1.5 GB

Equations

$$X_t:=X_0+{\int }_0^t{c}_{{\phi }_s}\mathrm{d}s;$$

$$Y_t:=Y_0+{\int }_0^t{r}_{{\phi }_s}(X_s)\mathrm{d}s.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{X}_t$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{Y}_t$$

$${\mathcal{F}}_i^{+}$$

$${\mathcal{F}}_i^{-}$$

$${\mathcal{F}}_i^0$$

$${\pi }_i(y)(\mathcal{A})$$

$$p_i(\mathcal{A})$$

$$\mathbf{\mu }=\left[\begin{array}{ccc}({\mu }_i^{+}{)}_{i\in {\mathcal{S}}_{+}}& ({\mu }_i^{-}{)}_{i\in {\mathcal{S}}_{-}}& ({\mu }_i^0{)}_{i\in {\mathcal{S}}_0}\end{array}\right],$$

$${\mu }_i^{\ell }(\mathcal{A}):={\mu }_i^{\ell }(\mathcal{A}\cap {\mathcal{F}}_i^{\ell })\text{mbox}for\mathcal{A}\subset \mathcal{F}.$$

$${\mu }_i^{\ell }{\mathbb{V}}_{ij}^{\ell m}(t)(\mathcal{A}):={\int }_{x\in {\mathcal{F}}_i^{\ell }}\mathrm{d}{\mu }_i^{\ell }(x)\text{mathds}P[{\phi }_t=j,X_t\in \mathcal{A}|{\phi }_0=i,X_0=x],$$

$$\mathbb{V}(t)=e^{\mathbb{B}t}\text{mbox}fort\ge 0.$$

$$\displaystyle\mathbb{B}=\left[\begin{array}[]{lll}\mathbb{B}^{++}&\mathbb{B}^{+-}&\mathbb{B}^{+0}\\ \mathbb{B}^{-+}&\mathbb{B}^{--}&\mathbb{B}^{-0}\\ \mathbb{B}^{0+}&\mathbb{B}^{0-}&\mathbb{B}^{00}\\ \end{array}\right],$$

$${\mu }_i^{\ell }(\mathcal{A})={\int }_{x\in \mathcal{A}}\mathrm{d}{\mu }_i^{\ell }(x)={\int }_{x\in \mathcal{A}}{\nu }_i^{\ell }(x)\mathrm{d}x+\text{mathds}{1}_{\left\{0\in \mathcal{A}\right\}}{p}_i^{\ell },$$

$$[\mathbf{\mu }\mathbb{B}(\mathcal{G}){]}_{\ell ,i}:=\sum _{m\in \{+,-,0\}}\sum _{j\in \mathcal{S}}{\mu }_i^{\ell }{\mathbb{B}}_{ij}^{\ell m}(\mathcal{G}\cap {\mathcal{F}}_j^m).$$

$$\displaystyle\mu_{i}^{\ell}\mathbb{B}_{ij}^{\ell m}(\mathcal{A})=\left\{\begin{array}[]{l}T_{ij}\left[\displaystyle\int_{x\in\mathcal{A}\cap\mathcal{F}^{\ell}_{i}}\nu_{i}^{\ell}(x)\,\mathrm{d}x+p_{i}^{\ell}\mathds{1}_{\left\{0\in\mathcal{A}\cap\mathcal{F}^{\ell}_{i}\right\}}\right]\quad\mbox{for }c_{j}\leq 0,\\ \vspace*{-0.2cm}\\ T_{ij}\left[\displaystyle\int_{x\in\mathcal{A}\cap\mathcal{F}^{\ell}_{i}}\nu_{i}^{\ell}(x)\,\mathrm{d}x+p_{i}^{\ell}\mathds{1}_{\left\{0\in\mathcal{F}^{\ell}_{i}\right\}}\mathds{1}_{\left\{u=0,\;v>0\right\}}\right]\quad\mbox{for }c_{j}>0.\end{array}\right.$$

$${\mu }_i^{\ell }{\mathbb{B}}_{ii}^{\ell m}(\mathcal{A})$$

$$\displaystyle=\left\{\begin{array}[]{l}-c_{i}\nu_{i}^{\ell}(v)\mathds{1}_{\{u\neq v\}}\mathds{1}_{\left\{v\in\partial_{L\backslash R}[\mathcal{F}^{\ell}_{i}]\right\}}-c_{i}\nu_{i}^{\ell}(0)\mathds{1}_{\{v=0\}}\mathds{1}_{\left\{0\in\partial_{L\backslash R}[\mathcal{F}^{\ell}_{i}]\right\}}\quad\mbox{for }c_{i}<0,\\ \vspace*{-0.2cm}\\ c_{i}\nu_{i}^{\ell}(u)\mathds{1}_{\{u\neq v\}}\mathds{1}_{\left\{u\in\partial_{R\backslash L}[\mathcal{F}^{\ell}_{i}]\right\}}\quad\mbox{for }c_{i}>0,\end{array}\right.$$

$${\mu }_i^{\ell }{\mathbb{B}}_{ii}^{\ell \ell }(\mathcal{A})$$

$$\displaystyle=\left\{\begin{array}[]{l}T_{ii}\left[\displaystyle\int_{x\in\mathcal{A}}\nu_{i}^{\ell}(x)\,\mathrm{d}x+p_{i}^{\ell}\mathds{1}_{\left\{0\in\mathcal{A}\right\}}\right]+\mathds{1}_{\{u\neq v\}}\left[c_{i}\nu_{i}^{\ell}(u)-c_{i}\nu_{i}^{\ell}(v)\mathds{1}_{\left\{v\not\in\partial_{R}[\mathcal{F}_{i}^{({\ell})}]\right\}}\right]\\ \vspace*{-0.2cm}\\ \hskip 28.45274pt-c_{i}\nu_{i}^{\ell}(0)\mathds{1}_{\{0\in\mathcal{A}\}}\mathds{1}_{\left\{0\not\in\partial_{R}[\mathcal{F}_{i}^{({\ell})}]\right\}}\quad\mbox{for }c_{i}<0,\\ \vspace*{-0.2cm}\\ T_{ii}\left[\displaystyle\int_{x\in\mathcal{A}}\nu_{i}^{\ell}(x)\,\mathrm{d}x\right]+\mathds{1}_{\{u\neq v\}}\left[c_{i}\nu_{i}^{\ell}(u)\mathds{1}_{\left\{u\not\in\partial_{L}[\mathcal{F}^{\ell}_{i}]\right\}}-c_{i}\nu_{i}^{\ell}(v)\right]\\ \vspace*{-0.2cm}\\ \hskip 28.45274pt\quad\mbox{for }c_{i}>0,\end{array}\right.$$

$${\mu }_i^{+}{\Psi }_{ij}(s)(\mathcal{A})$$

$$:={\int }_{x\in {\mathcal{F}}_i^{+}}\mathrm{d}{\mu }_i^{+}(x)\text{mathds}E\left[e^{-s\theta (0)}\text{mathds}{1}_{\left\{{\phi }_{\theta (0)}=j,X_{\theta (0)}\in \mathcal{A}\right\}}|X_0=x,Y_0=0,{\phi }_0=i\right].$$

$$\displaystyle\mathbb{U}(y,s)=\left[\begin{array}[]{cc}\mathbb{U}^{++}(y,s)&\mathbb{U}^{+-}(y,s)\\ \mathbb{U}^{-+}(y,s)&\mathbb{U}^{--}(y,s)\\ \end{array}\right],$$

$${\mu }_i^{\ell }{\mathbb{U}}_{ij}^{\ell m}(y,s)(\mathcal{A}):={\int }_{x\in {\mathcal{F}}_i^{\ell }}\mathrm{d}{\mu }_i^{\ell }(x)\text{mathds}E\left[e^{-s\omega (y)}\text{mathds}{1}_{\left\{{\phi }_{\omega (y)}=j,X_{\omega (y)}\in \mathcal{A}\right\}}|{\phi }_0=i,X_0=x\right].$$

$$\mathbb{U}(y,s)=e^{\mathbb{D}(s)y},$$

$${\mathbb{D}}_{ij}^{\ell m}(s)=[{\mathbb{R}}^{\ell }({\mathbb{B}}^{\ell m}-s\mathbb{I}+{\mathbb{B}}^{\ell 0}(s\mathbb{I}-{\mathbb{B}}^{00}{)}^{-1}{\mathbb{B}}^{0m}){]}_{ij},$$

$${\mu }_i^{\ell }{\mathbb{R}}_i^{\ell }(\mathcal{A}):={\int }_{x\in \mathcal{A}\cap {\mathcal{F}}_i^{\ell }}\frac{1}{r_i(x)}\mathrm{d}{\mu }_i^{\ell }(x).$$

$${\mathbb{D}}^{+-}(s)+\Psi (s){\mathbb{D}}^{-+}(s)\Psi (s)+{\mathbb{D}}^{++}(s)\Psi (s)+\Psi (s){\mathbb{D}}^{--}(s)=0.$$

$${\xi }_i(\mathcal{A}):=\underset{n\to \infty }{\lim }\text{mathds}P\left[X_{{\theta }_n}\in \mathcal{A},{\phi }_{{\theta }_n}=i\right].$$

$$\displaystyle\begin{bmatrix}\bm{\xi}&\bm{0}\end{bmatrix}\left(-\left[\begin{array}[]{ll}\mathbb{B}^{--}&\mathbb{B}^{-0}\\ \mathbb{B}^{0-}&\mathbb{B}^{00}\end{array}\right]\right)^{-1}\left[\begin{array}[]{l}\mathbb{B}^{-+}\\ \mathbb{B}^{0+}\end{array}\right]\Psi$$

$$\sum _{i\in {\mathcal{S}}_{-}}{\xi }_i({\mathcal{F}}_i^{-})$$

$$\mathbf{\pi }(y)=\left[\begin{array}{ccc}{\mathbf{\pi }}^{+}(y)& {\mathbf{\pi }}^{-}(y)& {\mathbf{\pi }}^0(y)\end{array}\right]=\left[\begin{array}{ccc}({\pi }_i^{+}(y){)}_{i\in {\mathcal{S}}_{+}}& ({\pi }^{-}(y{)}_{i\in {\mathcal{S}}_{-}}& ({\pi }^0(y){)}_{i\in {\mathcal{S}}_0}\end{array}\right],$$

$${\pi }_i^{\ell }(y)(\mathcal{A})={\pi }_i^{\ell }(y)(\mathcal{A}\cap {\mathcal{F}}_i^{\ell })\text{mbox}for\mathcal{A}\subset \mathcal{F}.$$

Full text

Here, we construct a discontinuous Galerkin method for approximating this stationary distribution, and illustrate the methodology using an on-off bandwidth sharing system, which is a special case of a stochastic fluid-fluid process.

1 Introduction

A stochastic fluid process $\{X_{t},\varphi_{t}\}_{t\geq 0}$ is a two-dimensional Markov process, where the phase $\varphi_{t}$ is a continuous-time Markov chain on a finite state space $\mathcal{S}$, and the fluid $X_{t}$ varies linearly at rate $c_{\varphi_{t}}$. A subset of Markov additive processes, stochastic fluids have been well-analysed in the past two decades. There have been two recent generalisations of stochastic fluid processes to a higher dimension: Miyazawa and Zwart [1] analysed discrete-time multidimensional Markov additive processes, and Bean and O’Reilly [2] studied the so-called stochastic fluid-fluid process, the latter is our focus in this paper.

A stochastic fluid-fluid is a Markov process $\{X_{t},Y_{t},\varphi_{t}\}_{t\geq 0}$, where the phase $\varphi_{t}$ is still a Markov chain on a finite state space $\mathcal{S}$; $X_{t}\in(-\infty,\infty)$ is the first fluid, which varies linearly at rate $c_{\varphi_{t}}$

[TABLE]

and $Y_{t}$ is the second fluid, which varies linearly at rate $r_{\varphi_{t}}(X_{t})$:

[TABLE]

As the classic fluid process $\{X_{t},\varphi_{t}\}_{t\geq 0}$ is used extensively in many areas, such as insurance and environmental modelling, it is clear that stochastic fluid-fluid models have even a wider range of applicability.

An example of application for a stochastic fluid fluid is the modelling of growth and bleaching of coral reefs, as described in [2]. In this process, we can model the density of symbiotic zooxanthellae at time $t$ by $X_{t}$, with the positive rates $c_{i}$ corresponding to the growth of the zooxanthellae, the negative rates to the bleaching. If the density $X_{t}$ is below a certain threshold $x$, the coral cannot store lipids until the density increases again past this level. During the time $X_{t}\in(0,x)$, the coral relies on stored lipids, modelled as $Y_{t}$, and dies when the latter runs out, that is, $Y_{t}=0$.

While the analyses in [1, 2] are markedly different, both papers drew inspiration from Neuts’ matrix-analytic approach [3, 4] to obtain the limiting behaviour of these processes, working with operators on function spaces instead of matrices. Thus, their closed-form expressions for the limiting distributions ([1, Theorem 4.1], [2, Theorem 2]) are given in terms of operators acting on measures, which are not immediately amenable to numerical computations for real-life applications. One way to numerically handle operators on function spaces is to construct approximations of the operators. To this end, there exist numerical procedures such as finite difference, finite volume, finite element, and discontinuous Galerkin (DG) methods [5]. The operators arising from fluid-fluid processes are assumed to be acting on a function space of smooth probability densities. The choice of an approximation method should reflect these properties in its solutions. In the DG method, the conservation of probability mass and local smoothness can be captured in the approximations [5].

In this paper, we construct a discontinuous Galerkin method to approximate the joint stationary distribution of a stochastic fluid-fluid process. We numerically illustrate the effectiveness of the methodology using an on-off bandwidth-sharing system of two processors [6]. In this example, inputs into the processors, $X_{t}$ and $Y_{t}$, are turned on and off by a Markov chain, $\varphi_{t}$; the combined output capacity is fixed and allocated according to the workload of the first, high-priority, processor $X_{t}$. Latouche et al. [6] evaluated the marginal limiting distribution of the first processor $X_{t}$, and provided bounds for the marginal limiting distribution of the workload of the second processor $Y_{t}$. We verify our DG approximations by comparing them against Monte Carlo simulations of the system, against analytical results obtained in [6], and against our intuitive understanding of the system dynamics. In all considered cases, we find the approximations to be accurate.

The paper is organised as follows. In Section 2, we give relevant background to present the joint stationary distribution of a stochastic fluid-fluid process. We construct in Section 3 a discontinuous Galerkin method to approximate the stationary distribution, and include numerical experiments in Section 4.

2 Preliminaries

Consider a stochastic fluid-fluid process $\{X_{t},Y_{t},\varphi_{t}\}_{t\geq 0}$. We assume that $X_{t},Y_{t}\in[0,\infty)$ and that there is a regulated boundary at level [math] for both buffers:

[TABLE]

for $i\in\mathcal{S}$. Let $T$ be the irreducible generator for the finite Markov chain $\varphi_{t}$. We denote by $C:=\mathrm{diag}(c_{i})_{i\in\mathcal{S}}$ the diagonal fluid-rate matrix for $X_{t}$, and $R(x):=\mathrm{diag}(r_{i}(x))_{i\in\mathcal{S}}$ the diagonal fluid-rate matrix for $Y_{t}$. For the remainder of this section, we summarize the findings of [2] on the joint stationary distribution of $\{X_{t},Y_{t},\varphi_{t}\}_{t\geq 0}$.

Let $\mathcal{F}$ be the state space of $X_{t}$, so $\mathcal{F}=[0,\infty)$. For each Markovian state $i\in\mathcal{S}$, we partition $\mathcal{F}$ according to the rates of change $r_{i}(\cdot)$ for the second fluid $Y_{t}$: $\mathcal{F}:=\mathcal{F}^{+}_{i}\cup\mathcal{F}^{-}_{i}\cup\mathcal{F}^{0}_{i},$ where

[TABLE]

For all $i\in\mathcal{S}$, the functions $r_{i}(\cdot)$ are assumed to be sufficiently well-behaved that $\mathcal{F}^{m}_{i}$, $m\in\{+,-,0\}$, is a finite union of intervals and isolated points. Moreover, define $\mathcal{S}_{+}:=\{i\in\mathcal{S}:\mathcal{F}^{+}_{i}\neq\varnothing\}$, $\mathcal{S}_{-}:=\{i\in\mathcal{S}:\mathcal{F}^{-}_{i}\neq\varnothing\}$, and $\mathcal{S}_{0}:=\{i\in\mathcal{S}:\mathcal{F}^{0}_{i}\neq\varnothing\}$.

We assume that the process $\{X_{t},Y_{t},\varphi_{t}\}$ is positive recurrent, in order to guarantee the existence of the joint stationary density operator $\bm{\pi}(y)=(\pi_{i}(y))_{i\in\mathcal{S}}$ and the joint stationary mass operator $\bm{p}=(p_{i})_{i\in\mathcal{S}}$, where for $\mathcal{A}\subset\mathcal{F}$

[TABLE]

The determination of $\bm{\pi}(y)$ involves two important matrices of operators, $\mathbb{B}$ and $\Psi$. Intuitively, for a set $\mathcal{A}\in\mathcal{F}$ and a measure vector $\bm{\mu}=(\mu_{i})_{i\in\mathcal{S}}$, $\bm{\mu}\mbox{e}^{\mathbb{B}t}(\mathcal{A})$ gives the conditional probability of $X_{t}\in\mathcal{A}$, and $\bm{\mu}\Psi(\mathcal{A})$ the conditional probability of $Y_{t}$ returning to level zero and doing so when $X_{t}\in\mathcal{A}$, given that the initial distribution is $\bm{\mu}$.

2.1 Matrix $\mathbb{B}$ of Operators

Let $\mathcal{M}(\mathcal{S}\times\mathds{R}_{+})$ be the set of integrable complex-valued Borel measures on the Borel $\sigma$-algebra $\mathcal{B}_{\mathcal{S}\times\mathds{R}_{+}}$. For $\bm{\mu}\in\mathcal{M}(\mathcal{S}\times\mathds{R}_{+})$, we can write

[TABLE]

where $\mu^{\ell}_{i}\in\mathcal{M}(\mathcal{F}^{\ell}_{i})$, the set of integrable complex-valued Borel measures on $\mathcal{B}_{\mathcal{F}_{i}^{\ell}}$, and

[TABLE]

We denote by $\mathbb{V}(t)$ the matrix of operators $\mathbb{V}_{ij}^{\ell m}(t):\mathcal{M}(\mathcal{F}_{i}^{\ell})\mapsto\mathcal{M}(\mathcal{F}^{m}_{j})$, $i\in\mathcal{S}_{\ell},j\in\mathcal{S}_{m}$ and $\ell,m\in\{+,-,0\}$, which are defined for $\mathcal{A}\subset\mathcal{F}_{j}^{m}$ as follows:

[TABLE]

the probability of the process $\{X_{t},\varphi_{t}\}$ being in the destination set $(\mathcal{A},j)$ at time $t$, given that it starts in $(\mathcal{F}^{\ell}_{i},i)$ according to the measure $\mu_{i}^{\ell}$. We can write $\mathbb{V}(t)$ in terms of its infinitesimal generator $\mathbb{B}$ as

[TABLE]

For $\mathbb{V}(t)$ and $\mathbb{B}$, as well as other operators of the same dimensions to be introduced later in the paper, the operators are partitioned according to $\{+,-,0\}$; for example,

[TABLE]

where each block $\mathbb{B}^{\ell m}$ is the $|\mathcal{S}_{\ell}|\times|\mathcal{S}_{m}|$ matrix of operators $\mathbb{B}^{\ell m}_{ij}$. We note that $\mathbb{V}(t)$ forms a strongly continuous semigroup on the set of measures that are absolutely continuous with respect to Lebesgue measure and have analytic densities as well as possibly a point mass at zero in certain phases [2]. Thus, we restrict the domain of $\mathbb{V}(t)$ to the set of such measures, and write

[TABLE]

where $\nu_{i}^{\ell}$ is the associated density and $p^{\ell}_{i}$ is the probability mass at the boundary [math] when $\varphi_{t}=i$, for $\ell\in\{+,-,0\}$ and $i\in\mathcal{S}$.

For simplicity, from here on we assume a set $\mathcal{A}$ is an interval, which might or might not include its end points, that is, $\mathcal{A}\in\{(u,v),(u,v],[u,v),[u,v]\}$. By [2, Lemma 3], the operators $\mathbb{B}_{ij}^{\ell m}:\mathcal{M}(\mathcal{F}^{\ell}_{i})\mapsto\mathcal{M}(\mathcal{F}^{m}_{j})$, $\ell,m\in\{+,-,0\}$ and $i,j\in\mathcal{S}$, are given as follows. We include also brief probabilistic interpretations of the terms (see remarks in [2] for more details). Note that in all of the following cases, we assume $\mathcal{A}\subset\mathcal{F}_{j}^{m}$; this is without loss of generality, as for any set $\mathcal{G}$ and a measure $\bm{\mu}$ in our chosen domain,

[TABLE]

1. Case 1.

When $i\neq j$,

[TABLE]

Case $1$ represents when there is a stochastic jump from state $i$ to state $j\neq i$, which happens with rate $T_{ij}$. The integral represents the probability mass of the intersection of the initiating domain $\mathcal{F}^{\ell}_{i}$ and the destination set $\mathcal{A}$. If $c_{j}\leq 0$, then the point mass $p_{i}^{\ell}$ is preserved after the change of phase. If $c_{j}>0$, then the point mass $p_{i}^{\ell}$ disperses into the density in state $j$, and is captured only if $\mathcal{A}$ has an upper bound strictly greater than [math]. Note that in this case [math] does not have to be in $\mathcal{A}$, it only has to be in the closure of $\mathcal{A}$.

1. Case 2.

When $i=j$ and $\ell\neq m$,

[TABLE]

where $\partial_{L\backslash R}[\mathcal{G}]$ denotes the left boundary point that mustn’t also be the right boundary point of the closure of the set $\mathcal{G}$, and similarly for $\partial_{R\backslash L}[\mathcal{G}]$.

Case $2$ represents a drift from $\mathcal{F}^{\ell}_{i}$ to $\mathcal{F}^{m}_{i}$. When neither $\mathcal{F}^{\ell}_{i}$ nor $\mathcal{A}$ is an isolated point, there is a transfer of density from one set to another through the relevant endpoints. When $\mathcal{A}=\{0\}$, $c_{i}<0$, and [math] is the left endpoint of the closure of $\mathcal{F}_{i}^{\ell}$ (which can’t be an isolated point itself), then there is also an accumulation of point mass.

1. Case 3.

When $i=j$ and $\ell=m$,

[TABLE]

where $\partial_{R}[\mathcal{G}]$ denotes the right boundary point of the closure of $\mathcal{G}$, and similarly for $\partial_{L}[\mathcal{G}]$. Case $3$ represents stochastic jumps out of state $i$ and drift across $\mathcal{F}^{\ell}_{i}$.

2.2 Matrix $\Psi$ of Operators

We denote by $\Psi(s)$ the $|\mathcal{S}_{+}|\times|\mathcal{S}_{-}|$ matrix of operators recording the Laplace-Stieltjes transforms of the time for $Y_{t}$ to return, for the first time, to the initial level of zero. Define the stopping time $\theta(y):=\inf\{t>0:Y_{t}=y\}$ to be the first time $Y_{t}$ hits level $y$, then each component $\Psi_{ij}(s):\mathcal{M}(\mathcal{F}^{+}_{i})\mapsto\mathcal{M}(\mathcal{F}^{-}_{j}),i\in\mathcal{S}_{+}$ and $j\in\mathcal{S}_{-}$, is given by

[TABLE]

Let $b(t):=\int_{0}^{t}\left|r_{\varphi_{z}}(X_{z})\right|\,\mathrm{d}z$ be the total unregulated amount of fluid that has flowed into or out of the second buffer $Y_{t}$ during $[0,t]$, and let $\omega(y):=\inf\{t>0:b(t)=y\}$ be the first time this accumulated in-out amount hits level $y$. We denote by $\mathbb{U}(y,s)$ the matrix of operators recording the Laplace-Stieltjes transforms of $\omega(y)$:

[TABLE]

where $\mathbb{U}^{\ell m}$ is the $|\mathcal{S}_{\ell}|\times|\mathcal{S}_{m}|$ matrix of operators $\mathbb{U}^{\ell m}_{ij}$, for $y>0$, $s\in\mathds{C}$, and Re$(s)\geq 0$. Each operator $\mathbb{U}_{ij}^{\ell m}(y,s):\mathcal{M}(\mathcal{F}_{i}^{\ell})\mapsto\mathcal{M}(\mathcal{F}^{m}_{j})$, for $\ell,m\in\{+,-\}$ and $i,j\in\mathcal{S}$, is given by

[TABLE]

We can write

[TABLE]

where $\mathbb{D}(s)$ is the infinitesimal generator of the strongly continuous semigroup $\mathbb{U}(\cdot,y)$. Lemma $4$ of [2] gives the following expression for $\mathbb{D}(s)$.

Lemma 2.1

For $y\geq 0$, $s\in\mathds{C}$ with Re$(s)\geq 0$, $\ell,m\in\{+,-\}$, and $i\in\mathcal{S}_{\ell},j\in\mathcal{S}_{m}$,

[TABLE]

where $\mathbb{R}^{\ell}:=\mathrm{diag}(\mathbb{R}_{i}^{\ell})_{i\in\mathcal{S}_{\ell}}$ is a diagonal matrix of operators $\mathbb{R}_{i}^{\ell}$ given by

[TABLE]

By [2, Theorem 1], $\Psi(s)$ has the following characterisation.

Theorem 2.2

For Re$(s)\geq 0$, $\Psi(s)$ satisfies the equation:

[TABLE]

Furthermore, if $s$ is real then $\Psi(s)$ is the minimal nonnegative solution.

2.3 Stationary Distribution

Let $\Psi:=\Psi(0)$. We define $\theta_{n}:=\inf\{t\geq\theta_{n-1}:Y_{t}=0\}$, for $n\geq 2$, to be the sequence of hitting times to level [math] of $Y_{t}$, with $\theta_{1}:=\theta(0)$. Consider a discrete-time Markov process $\{X_{\theta_{n}},\varphi_{\theta_{n}}\}_{n\geq 1}$, and for $i\in\mathcal{S}_{-}$ define the measure $\xi_{i}$ as follows

[TABLE]

By [2], the vector of measures $\bm{\xi}:=(\xi_{i})_{i\in\mathcal{S}}$ satisfies the following set of equations

[TABLE]

We reproduce Theorem 2 of [2] below, which gives the joint stationary distribution of $\{X_{t},Y_{t},\varphi_{t}\}$. Recall that the joint stationary density operator $\bm{\pi}(y)=(\pi_{i}(y))_{i\in\mathcal{S}}$ for $\{X_{t},Y_{t},\varphi_{t}\}$ and the joint stationary mass operator $\bm{p}=(p_{i})_{i\in\mathcal{S}}$ are defined by (4) and (5), respectively. We can write

[TABLE]

where

[TABLE]

Theorem 2.3

The density $\bm{\pi}^{\ell}(y)$, for $\ell\in\{+,-,0\}$ and $y>0$, and the probability mass $\bm{p}^{m}$, for $m\in\{-,0\}$, satisfy the following set of equations:

[TABLE]

where $\mathbb{K}:=\mathbb{D}^{++}(0)+\Psi\mathbb{D}^{(-+)}(0)$ and $\alpha$ is a normalizing constant.

3 Discontinuous Galerkin Approximations

Discontinuous Galerkin methods are used to approximate the solution to a system of partial differential equations. A brief description of these methods is as follows [5]. On the domain of the approximation, consider a finite sequence of so-called nodal points. We refer to each interval between two consecutive nodal points as a mesh, and the combination of meshes and nodal points as a stencil. Within each mesh, we have a finite element approximation, which constructs a finite-dimensional smooth Sobolev space by choosing appropriate piecewise polynomial basis functions, and then projects the partial differential equations onto this space. This projection leads to a new system of equations, referred to as the weak form of the original PDEs.

There is a flux operator moving probability from one mesh to another, in a manner similar to the underlying principle of a finite volume approximation: integrating the PDEs over each mesh and then constructing a new system of ordinary differential equations, which describe the change in the integral over the mesh. This method conserves probability, and can handle discontinuities, such as jumps and point masses.

Discontinuous Galerkin methods lead to global approximations in the space of piecewise functions. Intuitively, we sacrifice the continuity between meshes to gain the conservation of probability.

3.1 Application to a Stochastic Fluid-Fluid Model

Here, we construct a discontinuous Galerkin method to approximate the operator matrix $\mathbb{B}$ and subsequently the operator matrix $\Psi$, the two key ingredients of the joint stationary distribution for $\{X_{t},Y_{t},\varphi_{t}\}$. We begin with approximating the joint density $f_{i}(x,t)$ of $\{X_{t},\varphi_{t}\}$:

[TABLE]

which satisfies the system of partial differential equations

[TABLE]

subject to suitable boundary conditions [2].

While $X_{t}\in[0,\infty)$, any numerical approximation by necessity has to take place on a finite interval. Clearly, the state space truncation results in a point mass at the upper bound, which we have to address properly. It is important to choose an interval sufficiently large in order to control the error induced by the artificial upper bound for $X_{t}$. We shall further comment on this in Section 4, where we report our numerical experiments.

Let $[0,\mathcal{I}]$ be the domain of the approximation, where $\mathcal{I}<\infty$. We denote by $\left\{x_{1},x_{2},\ldots,x_{K}\right\}$ a finite sequence of $K$ nodal points on $[0,\mathcal{I}]$, with $x_{1}:=0$ and $x_{K}:=\mathcal{I}$, and by $\left\{\mathcal{D}_{1},\ldots,\mathcal{D}_{K-1}\right\}$ the sequence of corresponding meshes, $\mathcal{D}_{i}:=[x_{i},x_{i+1}]$, for $i=1,\ldots,K-1$. As there are two point masses, one at zero where there is a regulated boundary for $X_{t}$ and another at $\mathcal{I}$ where there is an artificial upper bound, the two end meshes, $\mathcal{D}_{1}$ and $\mathcal{D}_{K-1}$, each of which contains both a point mass and density.

For $k=1,\ldots,K-1$, we choose $N_{k}$ functions $\phi_{n}^{k}:\mathcal{D}_{k}\mapsto[0,\infty)$, $n=1,\ldots,N_{k}$, to be the basis functions, the span of which forms our approximation space, $V_{K}:=\oplus_{k=1}^{K-1}\{\phi_{1}^{k},\ldots,\phi_{N_{k}}^{k}\}$. Then, a function $u_{i}(\cdot,\cdot)\in V_{K}$ has the form:

[TABLE]

for some coefficient functions $\alpha^{k}_{i,n}(t)$. To construct an approximation for $f_{i}(x,t)$ we need to determine these functions $\alpha_{i,n}^{k}(t)$ or, equivalently, the $N$-dimensional row vector

[TABLE]

where

[TABLE]

and $N:=\sum\limits_{k=1}^{K-1}N_{k}$ is the total number of basis functions on $[0,\mathcal{I}]$.

To that end, there are three important matrices we need to introduce. The first two matrices, $M$ and $G$, are $N\times N$ block-diagonal and detail the dynamics within each mesh:

[TABLE]

where, for $m,n\in\{1,\ldots,N_{k}\}$,

[TABLE]

The third matrix, $F_{i}$, $i\in\mathcal{S}$, is related to the dynamics between adjacent meshes: it is the flux operator moving probability from one mesh to another. Let $u^{k}_{i}(x,t)$ be the projection of $u_{i}(x,t)$ onto the mesh $\mathcal{D}_{k}$. A central idea of the discontinuous Galerkin method is that the values of $u^{k}_{i}(x,t)$ on different meshes are linked to each other only through a numerical flux $f^{*}$ and that one can consider the approximation

[TABLE]

where $x_{k}^{R}$ and $x_{k}^{L}$ denote the right and left endpoints of the $k$th mesh, respectively, and $\left[g(x)\right]_{x_{a}}^{x_{b}}:=g(x_{b})-g(x_{a})$ for any function $g$.

There are many options for $f^{*}$. Here, we choose a first-order up-winding scheme [7], that is,

[TABLE]

for $i\in\mathcal{S}$ and $x\in\{x_{1}^{L},\ldots,x_{K}^{L},x_{1}^{R},\ldots,x_{K}^{R}\}$, where $\eta(\mbox{sgn}(c_{i}),x)$ is an adjustment parameter for when we have a stencil with meshes of different structures; note that the numerical flux $f^{*}_{i}(x,t)$ is defined at nodal points only. Suppose we are considering the $k$th mesh, $\mathcal{D}_{k}$. Then, we define the function $\eta(\mbox{sgn}(c_{i}),x)$ to be

[TABLE]

Here, the function $\eta_{\ell,k}$, $\ell\in\{k-1,k+1\}$, is the ratio of the integrals of the basis functions transferring the probability from the $\ell$th mesh to the $k$th mesh, where

[TABLE]

for one and then all $m\in\{1,\ldots,N_{k}\}$ and $n\in\{1,\ldots,N_{\ell}\}$. In other words, $\eta_{\ell,k}$ is not dependent on the particular choice of basis functions, $m$ and $n$, but on their meshes, $k$ and $\ell$, only.

This choice of $f^{*}$ requires the flux information only from the left if the fluid rate $c_{\varphi_{t}}$ of $X_{t}$ is positive, and only from the right if that rate is negative. Intuitively, when $c_{\varphi_{t}}>0$, the numerical flux of probability going from $\mathcal{D}_{k-1}$ into $\mathcal{D}_{k}$ is the probability density accumulated on the right-hand edge of the approximation on $\mathcal{D}_{k-1}$. Similarly, when $c_{\varphi_{t}}<0$, the numerical flux going from $\mathcal{D}_{k}$ into $\mathcal{D}_{k-1}$ is the probability density accumulated on the left-hand edge of the approximation on $\mathcal{D}_{k}$.

Thus, on nodal points we have

[TABLE]

where for any function $g$

[TABLE]

and is used to allow us access to the density on the left-hand edge of $\mathcal{D}_{k+1}$ and the right-hand edge of $\mathcal{D}_{k-1}$, respectively.

We are now ready to introduce the block-tridiagonal matrix $F_{i}$,

[TABLE]

where each sub-block $F_{i}^{\ell k}$ is of dimension $N_{\ell}\times N_{k}$ and thus

[TABLE]

is an $N$-dimensional row vector. Let $\bm{\phi}^{\ell}(x):=(\phi_{1}^{\ell}(x),\ldots,\phi^{\ell}_{N_{\ell}}(x))$ be the vector containing all basis functions on the $\ell$th mesh. We define for $c_{i}>0$

[TABLE]

Lemma 3.1

For $k=1,\ldots,K-1$ and $m=1,\ldots,N_{k}$,

[TABLE]

**Proof ** We begin with the RHS of (41). For $c_{i}>0$,

[TABLE]

As each basis function $\phi_{n}^{j}$, for $n=1,\ldots,N_{j}$, is trivially zero outside its $j$th mesh, we have

[TABLE]

Substituting (43) and (44) into (42) leads to

[TABLE]

The argument for $c_{i}<0$ follows analogously. \qed

Theorem 3.2

The weak formulation of the PDEs (26) is the following system of ordinary differential equations

[TABLE]

**Proof ** Consider Equation (26), for $i\in\mathcal{S},t\geq 0,x\in[0,\mathcal{I}]$, which we restate below:

[TABLE]

For details on the steps of discontinuous Galerkin methods, see [5]. Here, we start by replacing the density function $f_{i}(x,t)$ in (46) by its approximation $u_{i}(x,t)$ to obtain

[TABLE]

Multiplying both sides of (47) by a basis function $\phi_{m}^{k},m\in\{1,\ldots,N_{k}\}$ and $k\in\{1,\ldots,K=1\}$, and then integrating over the approximation domain $[0,\mathcal{I}]$ gives

[TABLE]

Expanding and then integrating the third term by parts leads to

[TABLE]

where the first part can be approximated by using the numerical flux $f^{*}_{i}(x,t)$, that is,

[TABLE]

Substituting (49) into (48) and expanding $u^{k}_{i}(x,t)$ gives

[TABLE]

Finally, we switch the order of integrals and summations to arrive at

[TABLE]

Equivalently, this can be written in matrix form as

[TABLE]

which completes the proof. \qed

Remark 3.3

Let $V_{K}$ be the approximation space as defined in (27). If

the weights $\bm{\alpha}_{i}(t)$ satisfy (45), 2. 2.

the eigenvalues of $c_{i}(G+F_{i})M^{-1}$ are in the negative real half of the complex plane including zero, and 3. 3.

only one basis function is non-zero on each boundary of each mesh,

then

[TABLE]

Intuitively, the second and third conditions are to ensure stability in the solution and conservation of the transfer of probability across meshes, respectively.

Defining the discontinuous Galerkin infinitesimal operator

[TABLE]

we construct a DG approximation $\mathcal{B}$ for the operator matrix $\mathbb{B}$ as follows. For $i\in\mathcal{S}$ and $\ell\in\{+,-,0\}$, define $\gamma_{i}^{\ell}$ to be an index set: given a region $\mathcal{F}_{i}^{\ell}$ of $X_{t}$ (defined in (1–3)), $\gamma_{i}^{\ell}$ is the set of meshes included in $\mathcal{F}^{\ell}_{i}$. For example, if $\mathcal{F}_{1}^{+}=\mathcal{D}_{1}\cup\mathcal{D}_{2}\cup\mathcal{D}_{5}$, then ${\gamma_{1}^{+}}=\{1,2,5\}$. Note that the nodal points should be chosen such that ${\gamma_{i}^{\ell}}\cap\gamma_{i}^{m}=\varnothing$ for all $\ell,m\in\{+,-,0\},\ell\neq m$, and $\gamma_{i}^{+}\cup\gamma_{i}^{-}\cup\gamma_{i}^{0}=\{1,2,\ldots,K-1\}$.

We choose each approximation matrix $\mathcal{B}_{ij}^{\ell m}$ to be $N\times N$. Similar to the definition for the operator matrix $\mathbb{B}$, introduced in Section 2.1, there are three cases.

1. Case 1.

When $i\neq j$, each $N_{k}\times N_{k}$ sub-block $\left[\mathcal{B}_{ij}^{\ell m}\right]_{kk}$ is given by

[TABLE] 2. Case 2.

When $i=j$ and $\ell\neq m$,

[TABLE] 3. Case 3.

When $i=j$ and $\ell=m$,

[TABLE]

Next, we approximate the operators $\mathbb{R}_{i}^{\ell}$ by an $N\times N$ block-diagonal matrix $\mathcal{R}_{i}^{\ell}$, with diagonal sub-blocks $[\mathcal{R}_{i}^{\ell}]_{kk}$ given by

[TABLE]

where $\rho_{i,(k,n)}\neq 0$ is an approximation of the fluid rate $r_{i}(X_{t})$ of $Y_{t}$, given $\varphi_{t}=i$, $X_{t}\in\mathcal{D}_{k}$, and the basis function $\phi_{n}^{k}$. These functions $\rho_{i,(k,n)}$ can sensibly be chosen to be:

[TABLE]

Define $\mathcal{B}^{\ell m}:=[\mathcal{B}_{ij}^{\ell m}]_{i\in\mathcal{S}_{\ell},j\in\mathcal{S}_{m}}$ for $\ell,m\in\{+,-,0\}$, and $\mathcal{R}^{\ell}:=\mathrm{diag}(\mathcal{R}^{\ell}_{i})_{i\in\mathcal{S}_{\ell}}$ for $\ell\in\{+,-\}$. Putting things together, a DG approximation of the operator $\mathbb{D}_{ij}^{\ell m}$ is given by the $N\times N$ matrix

[TABLE]

for $s\in\mathds{C}$ with Re$(s)>0$ and for $\ell,m\in\{+,-\}$.

Let $\mathcal{D}^{\ell m}(s):=[\mathcal{D}_{ij}^{\ell m}(s)]_{i\in\mathcal{S}_{\ell},j\in\mathcal{S}_{m}}$, then a DG approximation $\uppsi$ of the operator $\Psi$, the latter describing the probability of the fluid $Y_{t}$ starting at level zero and returning there for the first time, is an $(N\times|\mathcal{S}_{+}|)\times(N\times|\mathcal{S}_{-}|)$ matrix solution to the equation

[TABLE]

The matrix equation (52) can be solved using one of the many algorithms suggested by Bean et al. [8]. We then use this approximation $\uppsi$ to evaluate the limiting density, according to Theorem 2.3.

3.2 An Example

Consider a process $\{X_{t},Y_{t},\varphi_{t}\}$ where there is only one phase $i=1$ in $\mathcal{S}$, in other words, there is no Markov modulation. Consequently, the PDE (26) has the drift component only:

[TABLE]

Suppose we would like to develop a DG approximation for the density function $f(x,t)$ over an interval $[0,2.75]$. Consider the four meshes

[TABLE]

We choose the following basis functions:

[TABLE]

as depicted in Figure 1.

Then, we can verify that the matrices $M$ and $G$ are given by

[TABLE]

Assume that $c=1$. Then, the non-zero upper-diagonal blocks are given by

[TABLE]

with all other sub-blocks of $F$ being identically zero. Thus,

[TABLE]

where

[TABLE]

Consequently, the DG generator $\mathcal{Q}=c(G+F)M^{-1}$ is given by

[TABLE]

where note that all the row sums are zero.

Now, to illustrate the approximation $\mathcal{B}$ of the operator matrix $\mathbb{B}$, suppose the phase process $\varphi_{t}\in\mathcal{S}=\{1,2\}$ has a generator matrix $T$. Let the DG approximation for this fluid-fluid process $\{X_{t},Y_{t},\varphi_{t}\}$ have the same stencil as described in Figure 1.

When $\varphi_{t}=1$, we assume $Y_{t}$ has a positive fluid rate, $r_{1}(X_{t})>0$, when $X_{t}\in\mathcal{D}_{1}\cup\mathcal{D}_{2}$, and a negative rate, $r_{1}(X_{t})<0$, when $X_{t}\in\mathcal{D}_{3}\cup\mathcal{D}_{4}$. Thus, $\mathcal{F}^{+}_{1}=\mathcal{D}_{1}\cup\mathcal{D}_{2}$ and $\mathcal{F}^{-}_{1}=\mathcal{D}_{3}\cup\mathcal{D}_{4}.$ When $\varphi_{t}=2$, we assume the opposite: $\mathcal{F}^{-}_{2}=\mathcal{D}_{1}\cup\mathcal{D}_{2}$ and $\mathcal{F}^{+}_{2}=\mathcal{D}_{3}\cup\mathcal{D}_{4}$. Thus, $\gamma_{1}^{+}=\gamma_{2}^{-}=\{1,2\}$ and $\gamma_{1}^{-}=\gamma_{2}^{+}=\{3,4\}.$

1. Case 1.

When $i\neq j$: Suppose $i=1,j=2,\ell=+$, and $m=-$. Then,

[TABLE] 2. Case 2.

When $i=j$ and $\ell\neq m$: suppose $i=1,\ell=+$, and $m=-$. Then,

[TABLE]

If we are to pre-multiply $\mathcal{B}_{11}^{+-}$ with a row vector, we can see the matrix picks up the value of the right-most basis function with support in $\mathcal{F}_{1}^{+}$, and then distributes this to the first two basis functions with support in $\mathcal{F}_{1}^{-}$. This concurs with our physical interpretation (see Section 2.1) that there is a drift from $\mathcal{F}_{1}^{+}$ to $\mathcal{F}_{1}^{-}$. 3. Case 3.

When $i=j$ and $\ell=m$: Suppose $i=1$ and $\ell=-$. Then,

[TABLE]

4 Numerical Experiments

To illustrate the validity of our discontinuous Galerkin approximation, we perform numerical experiments on a stochastic fluid-fluid model, in a three-pronged approach.

First, we run Monte Carlo simulations, in order to compare the simulated joint density of $\{X_{t},\varphi_{t}\}$ evaluated at the time $Y_{t}$ first returns to the initial level [math] against that which is obtained via the return-probability matrix $\uppsi$. This numerically verifies the accuracy of our proposed approximation for the operator matrix $\Psi$. Second, using $\uppsi$ we evaluate the limiting joint density of $\{X_{t},\varphi_{t}\}$, which we compare against the same density analytically derived in [6]. Third, we vary the parameters of the second fluid $Y_{t}$ to confirm that the approximating joint density for $\{X_{t},\varphi_{t}\}$ does not change, while the marginal limiting distribution for $Y_{t}$ does, both of which are consistent with our intuitive understanding of the chosen example. In all three procedures, we find the approximations to be accurate.

We also analyse different choices for the level of spatial discretisation and the degree of polynomial basis functions, with respect to the order of convergence in relevant error terms.

4.1 An on-off bandwith-sharing model

The example we choose for our experiments is as follows. Consider a stochastic fluid-fluid $\{X_{t},Y_{t},\varphi_{t}\}_{t\geq 0},$ where $X_{t}$ and $Y_{t}$ represent the workloads in Buffer 1 and Buffer 2 at time $t\geq 0$, both driven by the phase $\varphi_{t},$ which is a Markov chain on the state space $\mathcal{S}=\{11,10,01,00\}$. Here, the state $11$ indicates inputs to both buffers being on, the state $00$ indicates both being off, the state $10$ is when only the first input is on, and the state $01$ is when only the second is on. The input of Buffer $k$ is switched from on to off with rate $\alpha_{k}$, and from off to on with rate $\beta_{k}$, for $k=1,2$. Thus, the infinitesimal generator $T$ for $\varphi_{t}$ is given by

[TABLE]

We denote by $\lambda_{k}$ the input rate of Buffer $k$ during an on period, and by $\zeta_{k}$ the output rate, and assume that $\lambda_{k}>\zeta_{k}$, for $k=1,2$. This example imitates the model considered in [6], where the two buffers share a fixed total output capacity $\kappa>0$.

In particular, the output rate $\zeta_{1}$ is constantly $\theta_{1}$, except for when the buffer is empty, that is, when $X_{t}=0$. On the other hand, $\zeta_{2}$ varies depending on the values of the buffers. More specifically, as Buffer 1 is considered high-priority, this buffer is allocated the entire output capacity $\kappa$ whenever exceeding a pre-determined threshold $x^{*}$, leaving $\zeta_{1}=\kappa$ and $\zeta_{2}=0$. However, when Buffer 1 is empty, Buffer 2 is given the whole output capacity $\kappa$, meaning $\zeta_{2}=\kappa$. For $X_{t}\in(0,x^{*})$, $\zeta_{1}=\theta_{1}$ and $\zeta_{2}=\theta_{2}$, for $\theta_{1},\theta_{2}\geq 0$ such that $\theta_{1}+\theta_{2}=\kappa$. Clearly, if Buffer $k$ is empty, its output rate $\zeta_{k}$ would be zero. Table 1 summarizes the output rates for all different scenarios.

Given its dynamics as currently defined based on [6], Buffer 1 is an example of what is known in the literature as a level-dependent fluid, because its net rates depend on the value of $X_{t}$ relative to the threshold $x^{*}$. As the existing theoretical analysis for stochastic fluid-fluid processes (developed in [2] and briefly summarized in Section 2) does not allow for level dependency in the first buffer, here we modify the bandwith-sharing model in [6] slightly. We let the output rate $\zeta_{1}$ of Buffer 1 remain $\theta_{1}$ for $X_{t}>0$, effectively eliminating the threshold effect on the first buffer but keeping the effect on the second. Table 2 represents the modified rates.

Consequently, the net rates of change for $X_{t}$, $c_{i}$, are given by

[TABLE]

and the net rates of change for $Y_{t}$, $r_{i}$, are as follows

[TABLE]

For our numerical experiments, we use the parameter choices given in [6]:

[TABLE]

As mentioned previously, while the true problem has an unbounded domain $[0,\infty)$, the discontinuous Garlekin method requires the domain of approximation to be a finite interval. Hence, for all approximations we consider a finite interval but large enough so that the boundary-induced dynamics do not significantly affect the results.

To specify the stencil for our numerical approximation, we define a vector \bm{\omega}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}K},h,\Delta_{h}} of $K$ nodal points as

[TABLE]

for $h>0$ and for $\Delta_{h}>0$, both sufficiently small. In this stencil, there are $K-1$ meshes, of which $K-5$ are interior meshes of length $h$, the left and right boundary meshes are of length $\Delta_{h}$, and the second-to-left and second-to-right ones are of length $h-\Delta_{h}$. The boundary meshes always have piecewise-constant approximations, because this is sufficient to approximate the point masses accumulated at either boundary.

4.2 The return-probability matrix $\uppsi$ via Monte Carlo simulations

We choose our initial distribution $v_{i}(x,y,0):=\mathds{P}\left[X_{0}\leq x,Y_{0}\leq y,\varphi_{0}=i\right]$ to be a point mass of $1$ for $(X_{0}=5,Y_{0}=0,\varphi_{0}=01)$, and zero everywhere else. The choice of an initial point being a point mass instead of a non-degenerate distribution is purely for the convenience of numerical simulation, because it eliminates the need of simulating multiple initial starting points. Using this initial condition and the parameters specified in (57, 58), we simulate $10^{5}$ trajectories, terminating each path either when Buffer $Y_{t}$ returns to zero or at time ${V}=10000$, Overall, 2.3% of the trajectories did not reach zero in buffer $Y$ by time $V$, and so are rejected.

Note that as we assume positive recurrence for our system, the probability of Buffer 2 returning to its initial level zero is $1$. However, the time this takes might be longer than the time window of our simulation, ${V}=10000$, which means the trajectories that we have terminated at time ${V}$ must eventually return to zero with probability $1$.

Let $\tau:=\inf\{t>0:Y_{t}=0\}$ be the first time $Y_{t}$ returns to zero. For each retained simulated trajectory, we record the values of $X_{\tau}$ and $\varphi_{\tau}$. For the states $10$ (the first input being on, the second being off) and $00$ (off-off), we present in Figure 2 the cumulative distributions of $X_{\tau}$, $\mathds{P}[X_{\tau}\leq x,Y_{\tau}=0,\varphi_{\tau}=i]$, with $i=00,10$, as determined by the simulations as well as by a piecewise linear DG approximation constructed from the approximating matrix $\uppsi$ of operator $\Psi$. In this piecewise linear DG approximation, we consider the approximation interval [0,\mathcal{I}]={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}[0,16]}, and the parameters of the stencil \bm{\omega}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}K},h,\Delta_{h}}, defined in (59), take the following values: $K=43,h=0.4$ and $\Delta_{h}=0.001$. The boundary meshes each have a constant basis function of value $1$; each $k$th interior mesh \mathcal{D}_{k}:={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}[x_{k},x_{k+1}]} has two piecewise linear basis functions

[TABLE]

for $x\in\mathcal{D}_{k}$ and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k=2,\ldots,K-2}. As can be seen in Figure 2, these two distributions are closely matched.

4.3 The marginal density of $\{X_{t},\varphi_{t}\}$

Since Buffer 1, $X_{t}$, is independent of Buffer 2, $Y_{t}$, we can use results from the existing literature on stochastic fluid flows to obtain the marginal limiting density $\bm{\chi}(x)=(\chi_{i}(x))_{i\in\mathcal{S}}$ of $X_{t}$:

[TABLE]

On the other hand, we can use the operator $\Psi$ to compute, based on Theorem 2.3, the joint limiting density $\pi_{i}(\mathcal{A},y)$, where

[TABLE]

Thus, we can approximate the joint limiting density $\pi_{i}(y)(\mathcal{A})$, via a discontinuous Galerkin approximation $\uppsi$ of $\Psi$, and consequently form an approximation of the marginal limiting density $\bm{\chi}(x)$. In particular, we evaluate an approximation of $\pi_{i}(y)([0,x))$, and then integrate the approximating function over the domain of $y$, and finally differentiate with respect to $x$; note that

[TABLE]

Let two vectors $\overline{\bm{\chi}}$ and $\widehat{\bm{\chi}}$ denote respectively the piecewise constant and piecewise linear DG approximations of $\bm{\chi}$, obtained via the corresponding approximations $\uppsi$. We use \bm{\omega}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}43},0.4,0.001} as our stencil for the DG approximation and the nodal points at which we evaluate the analytical density function $\bm{\chi}$. Define

[TABLE]

We present in Figure 3 the analytical density $\bm{\chi}$ at given nodal points, a piecewise constant DG approximation, and a piecewise linear DG approximation.

The approximations reconstruct the general shape of the density reasonably well: for both piecewise constant and piecewise linear, we see most of the probability being concentrated in the point mass at $X_{t}=0$, and then decaying as $x$ increases. We observe from the subplot inside the left plot of Figure 3 that the piecewise constant approximation $\overline{\chi}_{\textsc{off}}$ underestimates the point mass, and redistributes the difference over the rest of the state space of $X_{t}$. Hence, there is more mass in the tails of the densities. On the other hand, the piecewise linear approximation $\widehat{\bm{\chi}}$ appears to be very close to the analytical solution.

4.4 Sensitivity analysis of the dynamics of $Y_{t}$

To further confirm that the discontinuous Galerkin approximation $\uppsi$ of the operator $\Psi$ accurately captures the dynamics of $Y_{t}$, we vary the rates at which the input to this buffer switches on and off (denoted by $\alpha_{2}$ and $\beta_{2}$, respectively). As we modify these rates, we should see a change in the distribution of probability between

[TABLE]

the limiting marginal density of $X_{t}$ when $Y_{t}=0$, and

[TABLE]

that of $X_{t}$ when $Y_{t}>0$. On the other hand, the sum of these two densities should be identical in all the different meaningful scenarios of $\alpha_{2}$ and $\beta_{2}$; that is, the sum ${\chi}^{0}(x)+{\chi}^{+}(x)$ should remain fixed and be equal to the sum ${\chi}_{\textsc{on}}(x)+{\chi}_{\textsc{off}}(x)$, for all $x$.

To that end, we keep our stencil $\bm{\omega}_{43,0.4,0.001}$ and basis functions fixed, and compute the marginal limiting density for different values of $\alpha_{2}$ and $\beta_{2}$. The results coincide with what we expect from the dynamics in $Y_{t}$.

From Table 3, we observe that as $\alpha_{2}$ (the rate at which the input for $Y_{t}$ switches off) increases, so does the probability of $Y_{t}$ being empty. Furthermore, even though there are different amounts of probabilities in the two marginal densities, $\widehat{\chi}^{0}(x)$ and $\widehat{\chi}^{+}(x)$, their sum remains the same as the sum of the marginal limiting densities $\widehat{\chi}_{\textsc{on}}(x)$ and $\widehat{\chi}_{\textsc{off}}(x)$ for all calculated values of $x$ (data not shown here). These numerical results indicate that the dynamics of $Y_{t}$ are captured well by the DG approximations.

4.5 Errors of approximation

In a discontinuous Galerkin approximation, we choose the smoothness of the basis functions and the level of spatial discretisation. Once a particular selection is made, we project our operators into a finite-dimensional linear operator space corresponding to these choices (see Section 3). It has been shown that operators such as $\mathbb{B}$ (defined in Section 2.1) under a DG approximation have an error which converges at the order of $\mathcal{O}(h^{s}),$ where $h$ is the discretisation and $s$ is the degree of the basis [9].

However, this result cannot be easily translated across to the operator $\Psi.$ The DG approximation of the $\Psi$ operator is constructed by taking the DG approximation of the operators $\mathbb{B}$ and then solving the Riccati equation (52) using the approximate operators. Further, we then use this approximation of $\Psi$ to derive an approximation for the limiting density $\bm{\pi}$. With such a construction, it is not trivial to determine how the error propagates through the process of solving the Riccati equation, and then through further calculations to determine $\bm{\pi}$. Determining bounds for the approximation errors of $\Psi$ and $\bm{\pi}$, as functions of the discretisation and basis selection, is a topic for future research.

As a preliminary step in this direction, we empirically investigate how the approximation error of the marginal limiting density of the fluid $X_{t}$ (see Section 4.3) changes with respect to the choices of basis functions and the levels of discretisation. We begin by introducing our normed vector space in which we compare the different levels of discretisation and smoothness. Recall, from (59), that the left boundary mesh of our approximations is of length $\Delta_{h}$, and all but two interior meshes are of length $h$, with two compensating meshes of width $h-\Delta_{h}$. Let $[0,\mathcal{I}]$ be the interval on which we approximate our solution, then both the approximations and the analytical solution belong to the space $\mathcal{S}\times\mathcal{C}^{-1}([0,\mathcal{I}])$, where $\mathcal{C}^{-1}([0,\mathcal{I}])$ is the set of functions with countably many discontinuities. We choose the right boundary mesh to be a piecewise-constant function. Then, for any function $g:\mathcal{S}\times[0,\mathcal{I}]\mapsto\mathds{R}$, where $g(i,\cdot)\in\mathcal{C}^{-1}([0,\mathcal{I}])$ for $i\in\mathcal{S}$, we define the star seminorm as follows:

[TABLE]

Essentially, the star seminorm is an extension of the $L^{1}$ norm which incorporates our interpretation: that the left and right boundary meshes are treated as point masses. That is, we only study the total mass in the intervals $[0,\Delta_{h}]$ and $[\mathcal{I}-\Delta_{h},\mathcal{I}]$ and not the distribution over them.

We conduct two numerical experiments. The first is to understand the effects of choosing piecewise-linear basis functions over piecewise-constant, the second experiment is to understand the effects of treating a point mass as a density on a short interval.

In the first experiment, we choose $\Delta_{h}=10^{-6}$ and $\mathcal{I}=16$. We then consider the error between the approximation and the reference solution in the star seminorm. The notations $\widehat{\chi}_{\bm{\omega}}$ and $\overline{\chi}_{\bm{\omega}}$ denote, respectively, a piecewise-linear approximation with stencil $\bm{\omega}$ and a piecewise-constant approximation with stencil $\bm{\omega}$.

In the left plot in Figure 4, we observe that the piecewise-constant approximation has an error that scales approximately $\mathcal{O}({h^{0.88}}),$ with respective to the mesh size $h$, and the piecewise linear approximation has an error that scales approximately $\mathcal{O}({h^{1.84}})$. We observe that the coarsest piecewise-linear basis approximation has an error similar to that of the finest piecewise-constant approximation, which is two orders of magnitude finer.

In the second experiment, we fix $h=1.0,$ $\mathcal{I}=16$, and the basis functions to be piecewise-linear, and then we scale $\Delta_{h}$ and observe the trend in the star seminorm. The approximation error will plateau past the length where the error caused by $h$ and $\mathcal{I}$ dominate over the error gains of reducing the width of the boundary mesh. Hence, we subtract from this approximation error the approximation error of a finer approximation, $\widehat{\chi}_{19,1.0,0.005}$. The right plot in Figure 4 shows that the error scales approximately $\mathcal{O}({\Delta_{h}^{1.7}})$.

The numerical experiments above indicate that an increase in the degree of the basis functions would result in an increase in the order of convergence of the error in the star seminorm. However, this experiment is not taking into consideration that by increasing the degree of basis functions we need to increase the number of basis (elements) in the mesh. For example, the piecewise-constant has one basis element in each mesh, while the piecewise-linear has two. With respect to storage, the overall number of elements in the entire stencil is increasing linearly with the order of the elements.

In Table 4, we give some computational statistics of the approximations used in Figure 4 to aid in the understanding of the trade-offs between memory, computation time, mesh size, and error. We observe for a mesh size $h=0.5$ that the piecewise-constant approximation uses roughly half the number of elements and is twice as fast, compared to the piecewise linear approximation. However, the piecewise-linear approximation is two hundred times more accurate than the piecewise-constant approximation.

In our particular examples, if we were given a restriction to use no more than a prescribed amount of elements, then choosing a larger mesh size with piecewise-linear elements would be a better strategy than choosing a smaller mesh size and using piecewise-constant elements. The generalisation of these observations and the exploration of higher order basis functions is the focus for future research.

5 Conclusions

Finite Differences and Finite Volume methods have been used in the past to approximate operators that arise in stochastic processes. In principle, these methods approximate the operators by higher dimensional linear operators. In these methods, intuitive notions of mass conservation and positivity are captured; however, regularity is lost, making highly regular probability distributions computationally intensive.

We proposed the application of the discontinuous Galerkin method to approximate stochastic operators, with the intent that its ability to incorporate local regularity and maintain mass conservation, will lead to more accurate approximations and a reduction in computational effort. To demonstrate this, we applied the discontinuous Galerkin method to approximate all the operators needed to construct the joint stationary distribution of a stochastic fluid-fluid process.

The numerical results showed that the approximation of the stationary distribution arising from DG approximations of the operators is accurate and effective. We also verified that the operators and their dynamics were captured accurately. Furthermore, in our example, we observed that adding more regularity in the basis functions led to a significant decrease in computational effort. The DG method also enabled us to obtain other performance measures of stochastic fluid-fluid processes that are also analytically presented by operators.

Future work includes determining error bounds for the approximations of the operator $\Psi$ as well as of the stationary distribution in general, and a more thorough investigation of the computational effort as higher-order basis functions are used.

Acknowledgements

The authors acknowledge the financial support of the Australian Research Council (ARC) through the Discovery Grants DP110101663 and DP180103106. Bean, Nguyen, and O’Reilly also acknowledge the support of ACEMS (ARC Centre of Excellence for Mathematical and Statistical Frontiers).