A proton therapy model using discrete difference equations with an example of treating hepatocellular carcinoma

1.   Introduction

The convection-diffusion equation(CDE) governs the transmission of particles and energy caused by convection and diffusion. The CDE is given by

where $\omega$ denotes coefficient of viscosity and $\vartheta$ is the phase velocity respectively and both are positive. The Eq (1.1) is subject to the IC,

and the BCs,

where, $\phi$, $f_0$ and $f_1$ are given smooth functions.

Numerous computational procedures have been developed in the literature for the CDE. Chawla et al. ^[1] developed extended one step time-integration schemes for the CDE. Daig et al. ^[2] discussed least-squares finite element method for the advection-diffusion equation(ADE). Mittal and Jain ^[3] revisited the cubic B-splines collocation procedure for the numerical treatment of the CDE. The characteristics method with cubic interpolation for ADE was presented by Tsai et al. ^[4]. Sari et al. ^[5] worked on high order finite difference schemes for solving the ADE. Taylor-Galerkin B-spline finite element method for the one-dimensional ADE was developed by Kadalbajoo and Arora ^[6]. Kara and Zhang ^[7] presented ADI method for unsteady CDE. Feng and Tian ^[8] found numerical solution of CDE using the alternating group explicit methods with exponential-type. Dehghan ^[9] used weighted finite difference techniques for the CDE. Further, Dehghan ^[10] developed a technique for the numerical solution of the three-dimensional advection-diffusion equation.

A second-order space and time nodal method for CDE was conducted by Rizwan ^[11]. Mohebbi and Dehghan ^[12] presented a high order compact solution of the one-dimensional heat equation and ADE. Karahan ^[13,14] worked on unconditional stable explicit and implicit finite difference technique for ADE using spreadsheets. Salkuyeh ^[15] used finite difference approximation to solve CDE. Cao et. al ^[16] developed a fourth-order compact finite difference scheme for solving the CDEs. The generalized trapezoidal formula is used by Chawla and Al-Zanaidi ^[17] to solve CDE. Restrictive Taylor approximation has been used by Ismail et al. ^[18] to solve CDE. A boundary element method for anisotropic-diffusion convection-reaction equation in quadratically graded media of incompressible flow was studied by Salam et al. ^[19]. Azis ^[20] obtained standard-BEM solutions o two types of anisotropic-diffusion convection reaction equations with variable coefficients. The principle inspiration of this investigation is that the introduced scheme offers the solution as a piecewise adequately smooth continuous function enabling us to discover a numerical solution at any point in the solution domain. Moreover, it is simple to implement and has incredibly diminished computational expense. The refinement of the scheme with new approximations has elevated the accuracy of the scheme. The methodology used by von Neumann is utilized to affirm that the introduced scheme is unconditionally stable. The scheme is implemented to various test problems and the results are contrasted with the ones revealed in ^[25,26,27]. For further related studies, the interested reader is referred to ^[28,29] and references therein.

The remaining part of the paper is organized in the following sequence. The numerical scheme is presented in section 2 which is based on the cubic B-spline collocation method. Section 3 deals with the Scheme's stability and convergence analysis. The comparison of our numerical results with the ones presented in ^[25,26,27] is also presented in this section. The study's finding is summed up in section 4.

2.   Materials and method

Derivation of the scheme

Define $\Delta t = \frac{T}{N}$ to be the time and $h = \frac{d-c}{M}$ the space step sizes for positive integers $M$ and $N$. Let $t_n = n \Delta t, \; n = 0, 1, 2, ..., N$, and $z_j = jh, \; j = 0, 1, 2... M$. The solution domain $c\leq z\leq d$ is evenly divided by knots $z_{j}$ into $M$ subintervals $[z_{j}, z_{j+1}]$ of uniform length, where $c = z_{0} < z_{1} < ... < z_{n-1} < z_{M} = d$. The scheme for solving (1.1) assumes approximate solution $V(z, t)$ to the exact solution $v(z, t)$ to be ^[23]

where $D_j(t)$ are unknowns to be calculated and $B^3_j (z)$ ^[23] are cubic B-spline basis functions given by

Here, just $B^3_{j-1}(z), B^3_{j}(z)$ and $B^3_{j+1}(z)$ are last on account of local support of the cubic B-splines so that the approximation $v_j^n$ at the grid point $(z_j, t_n)$ at $n^{th}$ time level is given as

The time dependent unknowns $D_j^n(t)$ are determined using the given initial and boundary conditions and the collocation conditions on $B^3_j (z)$. Consequently, the approximations $v_j^n$ and its required derivatives are found to be

where ${\alpha_1 = \frac{1}{6}, } \; {\alpha_2 = \frac{4}{6}, } \; {\alpha_3 = \frac{1}{2h}, } \; {\alpha_4 = 0.}$

The new approximation ^[24] for $(v_j^n)_{zz}$ is given as

The problem (2.1) subject to the weighted $\theta$-scheme takes the form

where $h_j^{n} = \vartheta (v_j^n)_{zz}-\omega (v_j^n)_{z}$ and $n = 0, 1, 2, 3, ...$. Now utilizing the formula, $(v_j^n)_t = \frac{v_j^{n+1}-v_j^n}{k}$ in (2.6) and streamlining the terms, we obtain

Observe that $\theta = 0$, $\theta = \frac{1}{2}$ and $\theta = 1$ in the system (2.7) correspond to an explicit, Crank-Nicolson and a fully implicit schemes respectively. We use the Crank-Nicolson approach so that (2.7) is evolved as

Subtituting (2.4) and (2.5) in (2.8) at the knot $z_0$ returns

Substituting (2.4) and (2.5) in (2.8) produces

Subtituting (2.4) and (2.5) in (2.8) at the knot $z_M$ yields

From (2.9), (2.10) and (2.11), we acquire a system of $(M+1)$ equations in $(M+3)$ unknowns. To get a consistent system, two additional equations are obtained using the given boundary conditions. Consequently a system of dimension $(M +3)\times(M +3)$ is obtained which can be numerically solved using any numerical scheme based on Gaussian elimination.

Initial State: To begin iterative process, the initial vector $D^0$ is required which can be obtained using the initial condition and the derivatives of initial condition as follows:

The system (2.12) produces an $(M+3)\times(M+3)$ matrix system of the form

where,

$D^0 = [D_{-1}^0, D_{0}^0, D_{1}^0, ..., D_{M+1}^0]^T$ and $b = [\phi'(z_0), \phi(z_0), ..., \phi(z_M), \phi'(z_M)]^T$.

3.   Results

3.1.   Stability analysis

The von Neumann stability technique is applied in this section to explore the stability of the given scheme. Consider the Fourier mode, $D_j^n = \sigma^n e^{i \beta h j}$, where $\beta$ is the mode number, $h$ is the step size and $i = \sqrt{-1}$. Plug in the Fourier mode into equation (2.8) to obtain

where,

${\gamma_1 = \frac{k \vartheta}{24h^{2}}, }$

${\gamma_2 = \alpha_1- \frac{k \omega \alpha_3}{2} -\frac{k\vartheta}{3h^{2}}, }$

${\gamma_3 = \alpha_2+ \frac{k \omega\alpha_4}{2} +\frac{3k \vartheta}{4h^{2}}, }$

${\gamma_4 = \alpha_1+ \frac{k \omega \alpha_3}{2} -\frac{k \vartheta}{3h^{2}}, }$

${\gamma_5 = \alpha_1+ \frac{k \omega \alpha_3}{2} +\frac{k\vartheta}{3h^{2}}, }$

${\gamma_6 = \alpha_2- \frac{k \omega\alpha_4}{2} -\frac{3k \vartheta}{4h^{2}}, }$

${\gamma_7 = \alpha_1- \frac{k \omega \alpha_3}{2} +\frac{k \vartheta}{3h^{2}}.}$

Dividing equation (3.1) by $\sigma^n e^{i \beta j h}$ and rearranging, we obtain

Using $\cos(\beta h) = \frac{e^{i \beta h}+e^{-i \beta h}}{2}$ and $\sin(\beta h) = \frac{e^{i \beta h}-e^{-i \beta h}}{2i}$ in equation (3.2) and simplifying, we obtain

where, ${E_1 = \alpha_1 + \frac{k \vartheta}{3h^{2}}, } \; {F_1 = \frac{ k \omega \alpha_3}{2}, } \; {E_2 = \alpha_1-\frac{k \vartheta}{3h^{2}}, } \; {F_2 = \frac{ k \omega \alpha_3}{2}.}$

Notice that $\beta \in [-\pi, \pi]$. Without loss of generality, we can assume that $\beta = 0$, so that Eq (3.3) reduces to

which proves unconditional stability.

3.2.   Convergence analysis

In this section, we present the convergence analysis of the proposed scheme. For this purpose, we need to recall the following Theorem ^[21,22]:

Theorem 1. Let $v(z)\in C^4[c, d]$ and $c = z_{0} < z_{1} < ... < z_{n-1} < z_{M} = d$ be the partition of $[c, d]$ and $V^*(z)$ be the unique B-spline function that interpolates $v$. Then there exist constants $\lambda_i$ independent of $h$, such that

First, we assume the computed B-spline approximation to (2.1) as

To estimate the error, $\|v(z, t)-V(z, t)\|_\infty$ we must estimate the errors $\|v(z, t)-V^*(z, t)\|_\infty$ and $\|V^*(z, t)-V(z, t)\|_\infty$ separately. For this purpose, we rewrite the equation (2.8) as:

where, $v^* = v_j^{n+1}, (v^*)_{z} = (v_j^{n+1})_{z}, (v^*)_{zz} = (v_j^{n+1})_{zz}$ and $r(z) = v_j^{n}-\frac{1}{2}k \omega (v_j^n)_{z}+\frac{1}{2}k \vartheta (v_j^n)_{zz}$. Equation (3.4) can be written in matrix form as:

where, $R = ND^n+h$ and

,

$h = [f_0(t_{n+1}), 0, ...0, f_1(t_{n+1})]^T$, $D^n = [D_{-1}^n, D_{0}^n, D_{1}^n, ..., D_{M+1}^n]^T$,

${q_1 = \alpha_1- \frac{k \omega \alpha_3}{2} - \frac{7k\vartheta}{12h^{2}}, } \; {q_2 = \alpha_2+ \frac{k \omega \alpha_4}{2} +\frac{11k\vartheta}{8h^{2}}, } \; {q_3 = \alpha_1+ \frac{k \omega \alpha_3}{2} - \frac{7k\vartheta}{6h^{2}}, } \; {q_4 = \frac{7k\vartheta}{12h^{2}}, }$

${q_5 = \frac{k\vartheta}{4h^{2}}, } \; {q_6 = \frac{k\vartheta}{24h^{2}}, } \; {q_7 = \alpha_1+ \frac{k \omega \alpha_3}{2} +\frac{7k\vartheta}{12h^{2}}, } \; {q_8 = \alpha_2- \frac{k \omega \alpha_4}{2} -\frac{11k\vartheta}{8h^{2}}, }$

${q_9 = \alpha_1- \frac{k \omega \alpha_3}{2} +\frac{7k\vartheta}{6h^{2}}, } \; {q_{10} = \alpha_1- \frac{k \omega \alpha_3}{2} - \frac{7k \vartheta}{6h^{2}}, } \; {q_{11} = \alpha_1+ \frac{k \omega \alpha_3}{2} - \frac{7k \vartheta}{12h^{2}}, }$

${q_{12} = \alpha_1+ \frac{k \omega \alpha_3}{2} + \frac{7k \vartheta}{6h^{2}}, } \; {q_{13} = \alpha_1- \frac{k \omega \alpha_3}{2} +\frac{7 k \vartheta}{12h^{2}}.}$

If we replace $v^{*}$ by $V^{*}$ in (3.4), then the resulting equation in matrix form becomes

Subtracting (3.6) from (3.5), we obtain

Now using (3.4), we have

From (3.8) and Theorem (1), we have

where $M_1 = \lambda_0 h^2+\frac{k\omega}{2}\lambda_1 h+\frac{k\vartheta}{2} \lambda_2.$ It is obvious that the matrix $A$ is diagonally dominant and thus nonsingular, so that

Now using (3.9), we obtain

Let $a_{j, i}$ denote the entries of $A$ and $\eta_j, 0 \leq j \leq M+2$ is the summation of $jth$ row of the matrix $A$, then we have

From the theory of matrices we have,

where $a_{k, j}^{-1}$ are the elements of $A^{-1}$. Therefore

Substituting (3.13) into (3.11) we see that

where $M_2 = \frac{M_1}{|\xi_l|}$ is some finite constant.

Theorem 2. The cubic B-splines $\{B_{-1}, B_{0}, ....B_{M+1}\}$ defined in relation (2.2) satisfy the inequality

Proof. Consider,

Now for $z\in[z_{j+1}, z_{j+2}]$, we have

${|B^3_{j-2} (z)| \leq \frac{4}{6}}$

${|B^3_{j-1} (z)| \leq\frac{1}{6}, }$

${|B^3_{j} (z)| \leq \frac{4}{6}, }$

${|B^3_{j+1} (z)| \leq\frac{1}{6}.}$

Then, we have

as required.

Now, consider

Using (3.14) and Theorem 2, we obtain

Theorem 3. Let $v(z)$ be the exact solution and let $V(z)$ be the cubic collocation approximation to $v(z)$ then the provided scheme has second order convergence in space and

Proof. From Theorem 1, we have

From (3.16) and (3.17), we obtain

where $\mu = \lambda_0h^2+ \frac{5}{3} M_2.$

3.3.   Applications of numerical scheme

In this section, some numerical calculations are performed to test the accuracy of the offered scheme. In all examples, we use the following error norms

The numerical order of convergence $p$ is obtained by using the following formula:

where $L_\infty(n)$ and $L_\infty(2n)$ are the errors at number of partition $n$ and $2n$ respectively.

Example 1. Consider the CDE,

with IC,

and the BCs,

The analytic solution of the given problem is $v(z, t) = \exp(5z-(0.25+0.01\pi^2)t)\sin(\pi z)$. To acquire the numerical results, the offered scheme is applied to Example 1. The absolute errors are compared with those obtained in ^[25] at various time stages in Table 1. In Table 2, absolute errors and error norms are presented at time stages $t = 5, 10,100$. Figure 1 displays the comparison that exists between exact and numerical solutions at various time stages. Figure 2 depicts the 2D and 3D error profiles at $t = 1$. A 3D comparison between the exact and numerical solutions is presented to exhibit the exactness of the scheme in Figure 3.

The approximate solution when $t = 1, k = 0.01$ and $M = 20$ for Example 1 is given by

Example 2. Consider the CDE,

with IC,

and the BCs,

The analytic solution is $v(z, t) = \exp(0.22z-(0.0242+0.5\pi^2)t)\sin(\pi z)$. The numerical results are obtained by utilizing the proposed scheme. In Table 3, the comparative analysis of absolute errors is given with that of ^[25]. Absolute errors and errors norms at time levels $t = 5, 10,100$ are computed in Table 4. Figure 4 shows a very close comparison between the exact and numerical solutions at various stages of time. Figure 5 plots 2D and 3D absolute errors at $t = 1$. In Figure 6, a tremendous 3D contrast between the exact and approximate solutions is depicted.

The numerical solution when $t = 1, k = 0.01$ and $M = 20$ for Example 2 is given by

Example 3. Consider the CDE,

with IC,

and the BCs,

The analytic solution is $v(z, t) = \exp(0.25z-(0.0125+0.2\pi^2)t)\sin(\pi z)$. By utilizing the proposed scheme the numerical results are acquired. An excellent comparison between absolute errors computed by our scheme and the scheme of ^[25] is discussed in Table 5. In Table 6, absolute errors and error norms are computed at different time stages. A close comparison between the exact and numerical solutions at different time stages is depicted in Figure 7. Figure 8 plots 2D and 3D absolute errors at $t = 1$. Figure 9 deals with the 3D comparison that occurs between the exact and numerical solutions. The numerical solution when $t = 1, k = 0.01$ and $M = 20$ for Example 3 is given by

Example 4. Consider the CDE,

with IC,

and the BCs,

The analytic solution of the given problem is $v(z, t) = \sqrt{\frac{20}{20+t}}\exp(-\frac{(z-2-0.8t)^2}{0.4(20+t)})$. In Table 7, a comparison between absolute errors calculated by our scheme and the scheme of ^[26] is presented. Absolute errors and errors norms at time levels $t = 5, 10,100$ are presented in Table 8. Figure 10 illustrates the behavior of numerical solutions at various time stages. Figure 11 depicts the 2D and 3D graphs of absolute errors. In Figure 12, an excellent 3D contrast between the exact and numerical solutions shows the enormous accuracy of the scheme.

The approximate solution when $t = 1, k = 0.01$ and $M = 20$ for Example 4 is given by

Example 5. Consider the CDE,

with IC,

and the BCs,

The analytic solution of the given problem is $v(z, t) = \exp(1.17712434446770z-0.09t)$. The numerical outcomes are acquired by using the introduced scheme. The absolute errors and errors norms are computed in Tables 9 by applying the presented scheme on Example 5 and are compared with those obtained in ^[27]. In Table 10, absolute errors and error norms are presented at various time stages. Figure 13 plots the behavior of exact and numerical solutions at various time stages. Figure 14 depicts the 2D and 3D absolute error profiles at $t = 1$. In Figure 15, a 3D contrast between the exact and numerical solutions is presented and all the graphs are in good agreement.

The approximate solution when $t = 1, k = 0.01$ and $M = 20$ for Example 5 is given by

4.   Conclusions

This research uses a new approximation for second-order derivatives in the cubic B-spline collocation method to obtain the numerical solution of the CDE. The smooth piecewise cubic B-splines have been used to approximate derivatives in space whereas a usual finite difference has been used to discretize the time derivative. Special consideration is paid to the stability and convergence analysis of the scheme to ensure the errors do not amplify. The approximate solutions and error norms are contrasted with those reported previously in the literature. From this analysis, we can conclude that the estimated solutions are in perfect accord with the actual solutions. The scheme can be applied to a wide range of problems in science and engineering.

Acknowledgments

The authors are thankful to the worthy reviewers and editors for the useful and valuable suggestions for the improvement of this paper which led to a better presentation.

Conflict of interest

The authors have no conflict of interest.