Application of control theory in a delayed-infection and immune-evading oncolytic virotherapy

1.   Introduction

Telegraph equation is introduced by Oliver Heaviside and is a linear second-order hyperbolic partial differential equations that describe the current and voltage on an electrical transmission line with distance. The model demonstrates that wave patterns can form along the line and that the electromagnetic waves can be reflected on the wire. The nonhomogeneous telegraph equation with boundary and initial conditions is given by

with initial and Dirichlet boundary conditions

where $\alpha$ and $\beta$ are the real constants and $\Omega: = [0, 1]$. This equation referred to as the second-order hyperbolic partial differential equation with constant coefficients, models a mixture between wave propagation and diffusion by introducing a term that accounts for effects of finite velocity to the standard heat or mass transport equation. This equation represents a damped wave motion when $\alpha > 0$ and $\beta = 0$.

In this paper, we employ the multi-wavelet Galerkin method to solve nonhomogeneous telegraph equation (1.1) with initial (1.2) and boundary conditions (1.3). The Alpert's multi-wavelet bases are infinitely differentiable while having small compact support. In other words, these bases have combined the advantages of both finite difference and spectral bases. The first application of Alpert's multi-wavelet bases to the solution of PDEs is the adaptive solution of nonlinear time-dependent PDEs ^[1]. In this approach, the multiresolution representation of the derivative operator introduced, and then an adaptive solver developed for both linear and nonlinear PDEs. Further, multi-wavelet methods have been developed for PDEs such as conservation laws ^[10,16,23]. For other similar studies to related PDEs we refer to ^[4,5,6,7]

The required conditions for the existence of a unique solution to a nonhomogeneous telegraph equation with initial and boundary conditions an integral boundary condition via Galerkin's method investigated in ^[3,14]. The reproducing kernel Hilbert space method is utilised to solve this equation ^[3]. Lakestani et al. ^[19] employed a numerical solution based on the Galerkin and collocation method to solve this equation appropriately. In ^[18], the authors proposed the differential quadrature algorithm to obtain an approximate solution of the two-dimensional telegraph equation. A fast and simple method based on the Chebyshev wavelets method is proposed by Heydari et al., ^[15]. In this paper, the matrices of integration and differentiation are applied to reduce complexity. Mittal et al. ^[20] used cubic B-spline collocation method, whereas Dehghan and Shorki ^[11] proposed an algorithm based on thin plates spline radial basis functions using collocation points for solving this equation. A high accuracy method for the long-time evolution of the acoustic wave equation is introduced in ^[21]. Authors of ^[8] presented dual reciprocity boundary integral equation method. Due to the importance of this equation, many numerical methods have been proposed to solve the telegraph equation such as Quardatic B-spline collocation method ^[12], collocation method based on Chebyshev cardinal function ^[9], Lagrange interpolation and modified cubic B-spline differential quadrature methods ^[17], a hybrid meshless method ^[26], generalized finite difference method ^[25].

The paper is structured as follows. A brief introduction of the Alpert's multi-wavelets is provided in Section 2. In Section 3, the wavelet Galerkin method is used to approximate the solution of the problem, and the convergence analysis is investigated. Some numerical experiments are solved to illustrate the efficiency and accuracy of the proposed method in Section 4. finally conclusions are included in Section 5.

2.   Alpert's multi-wavelet

Assume that $\Omega: = [0, 1] = \cup_{b\in \mathcal{B}}X_{J, b}$ is the finite discretizations of $\Omega$, where $X_{J, b}: = [x_b, x_{b+1}], b\in\mathcal{B}: = \{0, \ldots, 2^J-1\}$ with $J\in \mathbb{Z}^+\cup \{0\}$, are determined by the point $x_b: = b/(2^J)$. On this discretization, appling the dilation $\mathcal{D}_{2^j}$ and the translation $\mathcal{T}_b$ operators to primal scaling functions $\{\phi_{0, 0}^0, \cdots, \phi_{0, 0}^{r-1}\}$, one can introduce the subspaces

of scaling functions. Here $\mathcal{R} = \{0, 1, \cdots, r-1\}$ and the primal scaling functions are the Lagrange polynomials of degree less than $r$ that introduced in ^[1].

Every function $p\in L^2(\Omega)$ can be represented in the form

where $\langle., . \rangle$ denotes the $L^2$-inner product

and $\mathcal{P}_J^r$ is the orthogonal projection that maps $L^2(\Omega)$ onto the subspace $V_J^r$. To find the coefficients $p_{J, b}^k$ that are determined by $\langle p, \phi_{J, b}^k\rangle = \int_{X_{J, b}}f(x)\phi_{J, b}^k(x)dx$, we shall compute these integrals. We apply the $r$-point Gauss-Legendre quadrature by a suitable choice of the weights $\omega_k$ and nodes ${\tau}_k$ for $k\in \mathcal{R}$ to avoid these integrals ^[1,24], via

Convergence analysis of the projection $\mathcal{P}_J^r(p)$ is investigated for the $r$-times continuously differentiable function $p\in \mathbb{C}^r(\Omega)$.

For the full proof of this approximation and further details, we refer the readers to ^[2]. Thus we can conclude that $\mathcal{P}_J^r(p)$ converges to $p$ with rate of convergence $O(2^{-Jr})$.

Let $\Phi_J^r$ be the vector function ${\Phi_J^r}: = [\Phi_{r, J, 0}, \cdots, \Phi_{r, J, 2^J-1}]^T$ and consists of vectors $\Phi_{r, J, b}: = [\phi_{J, b}^0, \cdots, \phi_{J, b}^{r-1}]$. The vector function $\Phi_J^r$ includes the scaling functions and called multi-scaling function. Furthermore, by definition of vector $P$ that includes entries $P_{br+k+1}: = p_{J, b}^k$, we can rewrite Eq (2.2) as follows

where $P$ is an $N$-dimensional vector $(N: = r2^J)$. The building blocks of these bases construction can be applied to approximate a higher-dimensional function. To this end, one can introduce the two-dimensional subspace $V_J^{r, 2}: = V_J^r\times V_J^r\subset L^2(\Omega)^2$ that is spanned by

Thus by this assumption, to derive an approximation of the function $p\in L^2(\Omega)^2$ by the projection operator $\mathcal{P}_J^r$, we have

where components of the square matrix ${P}$ of order $N$ are obtained by

where $\hat{\tau}_k = (\tau_k+1)/2$. Consider the $2r$-th partial derivatives of $f:\Omega^2\rightarrow \mathbb{R}$ are continuous. Utilizing this assumption, the error of this approximation can be bounded as follows

where $\mathcal{M}_{\max}$ is a constant.

By reviewing the spaces $V_J^r$, it is obvious these bases are nested. Hence there exist complement spaces $W_J^r$ such that

where $\oplus$ denotes orthogonal sums. These subspaces are spanned by the multi-wavelet basis

According to (2.8), the space $V_{J}$ may be inductively decomposed to $V_J^r = V_0^r\oplus (\oplus_{j = 0}^{J-1}W_j^r)$. This called multi-scale decomposition and spanned by the multi-wavelet bases and single-scale bases. This leads us to introduce the multi-scale projection operator $\mathcal{M}_J^r$. Assume that the projection operator $\mathcal{Q}_j^r$ the maps $L^2(\Omega)$ onto $W_j^r$. Thus we obtain

and consequently, any function $p\in L^2(\Omega)$ can be approximated as a linear combination of multi-wavelet bases

where

Note that, we can compute the coefficients $p_{0, 0}^k$ by using (2.2). But multi-wavelet coefficients from zero up to higher-level $J-1$ in many cases must be evaluated numerically. To avoid this problem, we use multi-wavelet transform matrix $T_J$, introduced in ^[22,24]. This matrix connects multi-wavelet bases and multi-scaling functions, via,

where ${{{\Psi}_J^r}}: = [\Phi_{r, 0, b}, \Psi_{r, 0, b}, \Psi_{r, 1, b}, \cdots, \Psi_{r, J-1, b}]^T$ is a vector with the same dimension $\Phi_J^r$ (here $\Psi_{r, j, b}: = [\psi_{j, b}^0, \cdots, \psi_{j, b}^{r-1}]$). This representation helps to rewrite Eq (2.10) as to form

where we have the $N$-dimensional vector ${\tilde{P}_J}$ whose entries are $p_{0, 0}^{k}$ and $\tilde{p}_{j, b}^k$ and is given by employing the multi-wavelet transform matrix $T_J$ as ${\tilde{P}_J} = T_J{P_J}$. Note that according to the properties of $T_J$ we have $T_J^{-1} = T^T_J$.

The multi-wavelet coefficients (details) become small when the underlying function is smooth (locally) with increasing refinement levels. If the multi-wavelet bases have $N_{\psi}^r$ vanishing moment, then details decay at the rate of $2^{-JN_{\psi}^r}$ ^[16]. Because vanishing moment of Alpert's multi-wavelet is equal to $r$, one can obtain ${\tilde{p}_{J, b}^k}\approx O(2^{-Jr})$ consequently. This allows us to truncate the full wavelet transforms while preserving most of the necessary data. Thus we can set to zero all details that satisfy a certain constraint $\varepsilon$ using thresholding operator $\mathcal{C}_{\varepsilon}$

and the elements of $\bar{P}_J$ are determined by

where $D_\varepsilon: = \{(j, b, k):|\tilde{p}_{j, b}^k| > \varepsilon\}$. Now we can bound the approximation error after thresholding via

where $\mathcal{P}_{J, D_{\varepsilon}}^r (p)$ is the projection operator after thresholding with the threshold $\varepsilon$ and $C_{thr} > 0$ is constant independent of $J, \varepsilon$.

3.   Multi-wavelet Galerkin method

Let us consider the generalized telegraph equations (TE) on the region $\Omega\times \Omega$ governed by the partial differential equation

with initial and Dirichlet boundary conditions

In order to derive the multi-wavelet Galerkin method for solving TE (3.1), we assume that the approximate solution $u$ can be expanded by the Alpert's multi-wavelet bases $\Psi_J^r$, i.e.,

Taking the first and second derivative with respect to $x$ and $t$ from both sides of the Eq (3.4), one can get

where the matrix $D_{\psi}$ is used to represents the derivative of multi-wavelet defined by ^[10,23].

Inserting (3.4) into (3.1) and employing (3.5) we obtain the residual via

where the vector $\tilde{F}$ is obtained the same way as the direction in (2.5) i.e.,

with $\tilde{F} = T_JFT_J^{-1}$. The Galerkin method requires $r_J^r$ to satisfy $\langle r_J^r, \psi_{j, b}^k\psi_{j', b'}^{k'}\rangle = 0$. Multiplying (3.6) by $\Psi_J^r(x)$ from left and ${\Psi_J^r}^T(t)$ from right and integrating, we end up with

where we employ orthonormality of multi-wavelet bases and the local support of these bases.

Equation (3.7) gives $(N^2)^2$ independent equations

We obtain $4N-4$ other equations from boundary conditions (3.1) and (3.2) via equations (3.3) and (3.5),

The problem becomes a system of linear equations with $N^2$ equations and $N^2$ unknowns,

where $\mathcal{U}$ and $\Upsilon$ are the vectorization of $U$ and $\tilde{F}$, respectively. It should be noted here that since most of the elements of the matrix $\mathcal{A}$ are zero, We use appropriate methods such as the generalized minimal residual method ($GMRES$). After solving this system the approximate solution is implicitly represented by (3.4).

Convergence analysis

To investigate the convergence analysis of the multi-wavelet Galerkin method, we put

subtracting this equation from (3.1), we obtain

where $e_J: = u-u_J$. Taking $L^2$-norm from both sides and using the triangle inequality yields

Now suppose that

where $\mathcal{E}$ is the $(N \times N)$ matrix and thus, one can write

where we utilize the orthonormality of multi-wavelet bases. According to the previous section, for any function $p$, when multi-wavelet bases hava high vanishing moments and the function $p$ is smooth, $\langle p, \Psi_J^r\rangle$ decays fast in $J\rightarrow \infty$. By means of vanishing moments of Alpert's multi-wavelets and the matrix norms inequalities, we get

Obviously, using (2.7), we can find

where $\eta: = \kappa\mathcal{M}\sqrt{N}$ with $\kappa: = \|D_{\psi}^2+2\alpha D_{\psi}+\beta^2 I_{N}+D_{\psi}^2\|_{\infty}$ and $\mathcal{M}$ is a constant. Consequently, $\|e\|_2\rightarrow 0$ when $J\rightarrow \infty$.

4.   Numerical results

To show the efficiency and accuracy of will employ the proposed method to obtain the approximate solution of the following examples. All of the computations have been done by Maple and MATLAB simultaneously.

Example 4.1. Assume the telegraph equation (3.1) with initial and boundary conditions (3.2) and (3.3). Let

and function

The exact solution of this equation is $u(x, t) = \tan(\frac{x+t}{2})$ ^[8,20].

In Table 1, with fixed $\alpha = 10$ and $\beta = 5$, choosing different refinement levels $J$ and multiplicity $r$ guarantees our convergence investigation. In Table 2, results are also compared with other methods ^[12,20] in terms of $L_{\infty}, L_2$ errors at different times. In ^[12,20], The space and time discretized by the rate of $0.001$ while for the proposed method this is $0.25$. In view of this and the results, our method is better much more than them. Taking $r = 7$ and $J = 2$, the approximate solution and $L_{\infty}$ errors are shown in Figure 1.

Example 4.2. In this example, we consider $u(x, t) = (x-x^2)t^2e^{-t}$ is the analytical solution of the telegraph equation (1.1) with initial and boundary conditions

and

$L_2, L_{\infty}$ errors are reported in Table 3 taking $r = 5$ ad $J = 3$. Results have been compared with the results of ^[11,20]. These results indicate that the proposed method solves this equation better than them. Time and space steps in these papers have been reported $\Delta t = 0.001$ and $h = 0.01$ while they are equal to $0.125$ in our simulation. The graph of numerical solution $u_J$ and $L_{\infty}$ error are shown in Figure 2 and the exact and approximate solution at different values of the time $t$ and space $x$ are plotted in Figure 3.

Example 4.3. Consider the Eq (1.1) with the right hand side function

and the initial and boundary conditions

The exact solution of this problem is given in ^[8,20], as

The effect of multiplicity parameter $r$ is show in Figure 4. Table 4, shows a comparison among the $L_{\infty}, L_2$ errors for the proposed method and other methods ^[12,20]. Given this table, we can find the proposed method very flexible and better than others. The effects of the refinement level $J$ and multiplicity parameter $r$ on $L_{\infty}, L_2$ errors are given in Table 5.

5.   Conclusions

In this study, the multi-wavelet Galerkin method was used to obtain an approximate solution of the telegraph equation. This method reduces the problem to a sparse system of linear equations, and then this system is solved by the GMRES method. The convergence analysis was investigated and some numerical tests were guaranteed it. Numerical experiments were shown the ability and flexibility of the proposed method in comparison to other methods.

Acknowledgments

This project was supported by Researchers Supporting Project number (RSP-2020/210), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

The writers state that they have no known personal relationships or competing financial interests that could have appeared to affect the work reported in this work.