Mathematical Analysis of Non-Isothermal Reaction–Diffusion Models Arising in Spherical Catalyst and Spherical Biocatalyst

Abstract

The theory of dynamical systems and their widespread applications involve, for example, the Lane–Emden-type equations which are known to arise in initial- and boundary-value problems with singularity at $\mathrm{the} \mathrm{time} t=0$. The main objective of this paper is to make use of some mathematical analytic tools and techniques in order to numerically solve some reaction–diffusion equations, which arise in spherical catalysts and spherical biocatalysts, by applying the Chebyshev spectral collocation method. The proposed scheme has good accuracy. The results are demonstrated by means of illustrative graphs and numerical tables. The accuracy of the proposed method is verified by a comparison with the results which are derived by using analytical methods.

1. Introduction

The Lane–Emden-type equations are known to arise in initial- and boundary-value problems with singularity at $\mathrm{the} \mathrm{time} t=0$. Analytical solutions in the neighbourhood of $t=0$ are always possible to find (see [1,2]), such as those given in [3,4,5,6]:

$$v^{″}\left(t\right)+\frac{\gamma }{t}{v}^{\prime }\left(t\right)+g\left(v\left(t\right)\right)=0, 0\le t\le 1,\gamma >0$$

(1)

with the following boundary conditions:

$$v^{\prime }\left(0\right)=b_1 \mathrm{and} v\left(1\right)=b_2$$

(2)

Here, $v\left(t\right)$ denotes the unknown function on $\left[0,1\right]$, $v^{\prime }\left(0\right)$ denotes the derivative of $v\left(t\right)$ at time $t=0$, and $b_1 and b_2$ are constants.

Equation (1) arises in astrophysics and many other branches of science such as physics, chemistry and bio-mathematics (see, for example, [7,8,9,10,11]).

1.1. Spherical Catalyst Model

The Lane–Emden boundary value problem (BVP) models the chemical species dimensionless concentration within a spherical catalyst and is given in [12] as follows:

$$v^{″}\left(t\right)+\frac{2}{t}{v}^{\prime }\left(t\right)-{\rho }^{2}v\left(t\right)\exp \left(\frac{\alpha \mu \left(1-v\left(t\right)\right)}{1+\mu \left(1-v\left(t\right)\right)}\right)=0$$

(3)

with the following boundary conditions:

$$v^{\prime }\left(0\right)=0 \mathrm{and} v\left(1\right)=1$$

(4)

The effectiveness factor τ is defined by:

$$\tau =\frac{3}{{\rho }^2}{\frac{dv}{dt}|}_{t=1}$$

1.2. Mathematical Model of the Spherical Biocatalyst Equation

The chemical species dimensionless concentration within a spherical biocatalyst is modelled by the Lane–Emden equation and is written as follows (see [13]):

$$v^{″}\left(t\right)+\frac{2}{t}{v}^{\prime }\left(t\right)-{\rho }^2\frac{\left(1+\mu \right)v\left(t\right)}{1+\mu v\left(t\right)}=0$$

(5)

together with following boundary conditions:

$$v^{\prime }\left(0\right)=0 \mathrm{and} v\left(1\right)=1$$

(6)

Here, in this article, we consider spherical catalytic and spherical biocatalyst pellets, in which non-isothermal single reactions take place.

In refs. [14,15], analytical techniques are used to solve these models. In the year 2017, Wazwaz [16] used a variational iteration method to solve these models. Subsequently, in the year 2018, the OHAM technique was used to solve these spherical models (see, for details, [17]). By using the third-order approximation, these models were solved numerically in [18]. Furthermore, in [14], these spherical models are solved by using a numerical technique.

Here, in our present investigation, we propose a spectral collocation method by using the shifted Chebyshev polynomials for finding the approximate solutions of the spherical catalyst and spherical biocatalyst models. In this method, we approximate the unknowns by using the shifted Chebyshev polynomials and truncate then to a finite dimension. Then, by the use of approximations, the spherical catalyst and spherical biocatalyst models are changed into a system of simpler equations whose solution leads to the approximate solution of the spherical catalyst and spherical biocatalyst models. The spectral collocation method has been used by many authors in order to analyse mathematical models, including fractional-order models, with applications in physical and engineering sciences and in many other fields (see [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38]). Detailed applications of fractional calculus in science and engineering can be found in [39]. Some applications of iterative and spectral methods to solve real-world problems can be seen in [40,41,42,43]. The results, which we have found in this article, are demonstrated numerically and illustrated graphically. The main finding of this paper is that the proposed technique is easy for computation purposes and timesaving. The results will be helpful for researchers working on applied chemistry. We have performed all computations on the programming software MATLAB R2018b.

2. Preliminaries

The third-kind Chebyshev polynomials are defined by the following equivalent forms [44,45]:

$$V_m\left(t\right)=\frac{\cos \left(m+\frac{1}{2}\right)\theta }{\cos \frac{\theta }{2}}=\frac{2^{2m}}{\left(\begin{array}{c}2m\\ 2\end{array}\right)}{P}_m^{\left(-\frac{1}{2},\frac{1}{2}\right)}\left(t\right)$$

(7)

where $t=\cos \theta , \theta \in \left[0,\pi \right]$ and $P_m^{\left(\alpha , \beta \right)}\left(t\right)$ is the Jacobi polynomial of degree $m$.

The equivalent analytical form of shifted Chebyshev polynomials of the third kind ($\alpha =-\frac{1}{2}, \beta =\frac{1}{2}$) is given by

$$V_i\left(t\right)={\displaystyle \sum }_{k=0}^i{\left(-1\right)}^{i-k}\frac{\mathsf{\Gamma }\left(\mathrm{i}+\frac{3}{2}\right)\mathsf{\Gamma }\left(\mathrm{i}+\mathrm{k}+1\right)}{\mathsf{\Gamma }\left(\mathrm{k}+\frac{3}{2}\right)\mathsf{\Gamma }\left(\mathrm{i}+1\right)\left(\mathrm{i}-\mathrm{k}\right)!\mathrm{k}!}{t}^k$$

(8)

The shifted Chebyshev polynomials of the third kind are orthogonal in the interval [0, 1] with respect to the following weight function:

$$w\left(t\right)=\sqrt{\frac{t}{1-t}}$$

and have their orthogonality property as follows:

$${{\displaystyle \int }}_0^1{V}_n\left(t\right)V_m\left(t\right)w\left(t\right)dt=\{\begin{array}{c}\frac{\pi }{2}, n=m\\ 0, n\ne m\end{array}$$

(9)

A function $f\in L_{w\left(t\right)}^2\left[0, 1\right]$, with $\left|f^{″}\left(t\right)\right|\le Q$, can be written as follows:

$$f\left(t\right)=\underset{n\to \infty }{\lim }{\displaystyle \sum }_{i=0}^n{c}_i{V}_i\left(t\right),$$

(10)

where

$$c_i=\langle f\left(t\right),V_i\left(t\right)\rangle =\frac{2}{\pi }w\underset{0}{\overset{1}{{\displaystyle \int }}}f\left(t\right)V_i\left(t\right)w\left(t\right)dt$$

(11)

The finite-dimensional approximation for Equation (6) is given by

$$f\cong {\displaystyle \sum }_{i=0}^m{c}_i{V}_i\left(t\right)=C^T{\lambda }_m\left(t\right),$$

(12)

where

$$C={\left[c_0,c_1,\dots .,c_m\right]}^T\mathrm{and} {\lambda }_m\left(t\right)={\left[V_0,V_1,\dots .,V_m\right]}^T$$

(13)

Theorem 1.

If${\lambda }_m\left(t\right)={\left[V_0,V_1,\dots .,V_m\right]}^T$denotes the shifted Chebyshev vector and if$v>0,$then:

$$I^v{V}_i\left(t\right)=I^{\left(v\right)}{\lambda }_m\left(t\right),$$

(14)

where$I^{\left(v\right)}=\left(\epsilon \left(i,j\right)\right)$is the$\left(m+1\right)\times \left(m+1\right)$operational matrix of integral of order$v$and its$\left(i,j\right)$th entry is given by

$$\begin{array}{l}\epsilon \left(i,j\right)\\ ={\displaystyle \sum }_{k=0}^i{\displaystyle \sum }_{l=0}^j{\left(-1\right)}^{i+j-k-l}\frac{\mathsf{\Gamma }\left(\frac{1}{2}\right)\mathsf{\Gamma }\left(i+\frac{3}{2}\right)\mathsf{\Gamma }\left(i+k+1\right)\mathsf{\Gamma }\left(j+l+1\right)\mathsf{\Gamma }\left(v+k+l+\frac{3}{2}\right)\left(2j+1\right)j!}{\left(i-k\right)!\left(j-l\right)!\left(l\right)! \mathsf{\Gamma }\left(k+\frac{3}{2}\right)\mathsf{\Gamma }\left(i+1\right)\mathsf{\Gamma }\left(v+k+1\right)\mathsf{\Gamma }\left(j+\frac{1}{2}\right)\mathsf{\Gamma }\left(l+\frac{3}{2}\right)\mathsf{\Gamma }\left(k+l+v+2\right)}.\end{array}$$

(15)

Proof. 

Please see [28].

In this paper, we will use this matrix for $v=1, 2$. □

3. Outline of the Proposed and Used Method

In this section, we give an outline of our proposed method. Models with undetermined constants as an intermediate step in the determination of the non-linear BVP are written as follows:

$$v^{″}\left(t\right)+\frac{2}{t}{v}^{\prime }\left(t\right)-{\rho }^{2}v\left(t\right)\exp \left(\frac{\alpha \mu \left(1-v\left(t\right)\right)}{1+\mu \left(1-v\left(t\right)\right)}\right)=0$$

(16)

$$v^{″}\left(t\right)+\frac{2}{t}{v}^{\prime }\left(t\right)-{\rho }^2\frac{\left(1+\mu \right)v\left(t\right)}{1+\mu v\left(t\right)}=0$$

(17)

with

$$v\left(0\right)=\beta , v^{\prime }\left(0\right)=0 \mathrm{and} v\left(1\right)=1$$

(18)

where $v\left(0\right)=\beta$ represents an undetermined constant.

We first choose the following approximation:

$$v^{″}\left(t\right)=C^T{\lambda }_n\left(t\right)$$

(19)

Then, by the first-order integration of Equation (18), we obtain

$$v^{\prime }\left(t\right)=C^T{I}^{\left(1\right)}{\lambda }_n\left(t\right)+P^T{\lambda }_n\left(t\right)$$

(20)

Taking the first-order integration of Equation (19), we find that

$$v\left(t\right)=C^T{I}^{\left(2\right)}{\lambda }_n\left(t\right)+P^T{I}^{\left(1\right)}{\lambda }_n\left(t\right)+Q^T{\lambda }_n\left(t\right)$$

(21)

where

$$v^{\prime }\left(0\right)=P^T{\lambda }_n\left(t\right) \mathrm{and} v\left(0\right)=Q^T{\lambda }_n\left(t\right)$$

(22)

We now take

$$1=D^T{\lambda }_n\left(t\right)$$

(23)

3.1. The Spherical Catalyst Equation

Combining Equations (15) and (17)–(22), we obtain

$$C^T{\lambda }_n\left(t\right)+\frac{2}{t}\left(C^T{I}^{\left(1\right)}+P^T\right){\lambda }_n\left(t\right)-{\rho }^2\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(t\right)\exp \left(\frac{\alpha \mu \left(D^T-\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right)\right){\lambda }_n\left(t\right)}{(D^T+\mu (D^T-(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T))){\lambda }_n\left(t\right)}\right)=0$$

(24)

For Equation (23), the residual is defined by

$$R_n\left(t\right)=C^T{\lambda }_n\left(t\right)+\frac{2}{t}\left(C^T{I}^{\left(1\right)}+P^T\right){\lambda }_n\left(t\right)-{\rho }^2\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(t\right)\exp \left(\frac{\alpha \mu \left(D^T-\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right)\right){\lambda }_n\left(t\right)}{(D^T+\mu (D^T-(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T))){\lambda }_n\left(t\right)}\right)$$

(25)

Collocating Equation (24) at $n-1$ points, which are given by

$$t_i=\frac{i}{n}, i=1,2,\dots ,n-1$$

we have

$$R_n\left(t_i\right)=C^T{\lambda }_n\left(t_i\right)+\frac{2}{t_i}(C^T{I}^{\left(1\right)}+P^T){\lambda }_n\left(t_i\right)-{\rho }^2(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T){\lambda }_n\left(t_i\right)\exp (\frac{\alpha \mu \left(D^T-\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right)\right){\lambda }_n\left(t_i\right)}{(D^T+\mu (D^T-(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T))){\lambda }_n\left(t_i\right)})$$

(26)

The corresponding initial and boundary conditions can be written as follows:

$$\begin{gathered} \left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(0\right)=\beta \\ \left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(1\right)=1 \end{gathered}$$

(27)

From Equations (25)–(27), we obtain a system of $n+1$ equations. Solving this system, we obtain our solution for the spherical catalyst model.

3.2. The Spherical Biocatalyst Equation

Combining Equations (16)–(22), we obtain

$$C^T{\lambda }_n\left(t\right)+\frac{2}{t}\left(C^T{I}^{\left(1\right)}+P^T\right){\lambda }_n\left(t\right)-{\rho }^2\frac{\left(1+\mu \right)\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(t\right)}{\left(D^T+\mu \left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right)\right){\lambda }_n\left(t\right)}=0$$

(28)

For Equation (28), the residual is defined by

$$R_n\left(t\right)=C^T{\lambda }_n\left(t\right)+\frac{2}{t}\left(C^T{I}^{\left(1\right)}+P^T\right){\lambda }_n\left(t\right)-{\rho }^2\frac{\left(1+\mu \right)\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(t\right)}{\left(D^T+\mu \left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right)\right){\lambda }_n\left(t\right)}$$

(29)

Collocating Equation (28) at $n-1$ points which are given by

$$t_i=\frac{i}{n}, i=1,2,\dots ,n-1$$

we find that

$$\begin{array}{l}{R}_n\left(t_i\right)=C^T{\lambda }_n\left(t_i\right)+\frac{2}{t_i}\left(C^T{I}^{\left(1\right)}+P^T\right){\lambda }_n\left(t_i\right)\\ -{\rho }^2\frac{\left(1+\mu \right)\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(t_i\right)}{\left(D^T+\mu \left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right)\right){\lambda }_n\left(t_i\right)}\end{array}$$

(30)

The corresponding initial and boundary conditions can be written as follows:

$$\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(0\right)=\beta$$

(31)

$$(\left(C^T{I}^{\left(2\right)}+P^T{I}^{\left(1\right)}+Q^T\right){\lambda }_n\left(1\right)=1$$

(32)

From Equations (30)–(32), we obtain a system of $n+1$ equations. Solving this system, we obtain our solution for the spherical biocatalyst model.

4. Convergence Analysis

Theorem 2.

Let$\frac{d^{2}v\left(t\right)}{d{t}^2}\in L_{w\left(t\right)}^2\left[0,1\right]$and${\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n$be the$n^{th}$approximation of it obtained by$\left(n+1\right)$shifted Chebyshev polynomials of the third kind. Additionally, let$\left|\frac{d^{4}v\left(t\right)}{d{t}^4}\right|<M$for a positive constant M.Then, as$n\to \infty$, the approximated value${\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n$converges to the exact value$\frac{d^{2}v\left(t\right)}{d{t}^2}$with the following inequality for the expansion coefficients of$\frac{d^{2}v\left(t\right)}{d{t}^2}$ :

$$\left|c_n\right|<\frac{M}{n^2}$$

(33)

Proof. 

Let

$$\frac{d^{2}v\left(t\right)}{d{t}^2}={\displaystyle \sum }_{i=1}^{\infty }{c}_i{V}_i\left(t\right)$$

(34)

Taking the $n^{th}$ approximation of the above equation, we obtain

$${\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n={\displaystyle \sum }_{i=1}^n{c}_i{V}_i\left(t\right)$$

(35)

By subtracting Equation (35) from Equation (34), we find that

$$\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)-{\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n={\displaystyle \sum }_{i=n+1}^{\infty }{c}_i{V}_i\left(t\right)$$

(36)

The coefficients in Equation (36) are given by

$$\begin{gathered} c_n=\frac{2}{\pi }\underset{0}{\overset{1}{{\displaystyle \int }}}\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)V_n\left(t\right)w\left(t\right)dt \\ {c}_n=\frac{2}{\pi }\underset{0}{\overset{1}{{\displaystyle \int }}}\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)V_n\left(t\right)\sqrt{\frac{t}{1-t}}dt \end{gathered}$$

(37)

If we set $2t-1=\cos y$, then we have

$$c_n=\frac{2}{\pi }\underset{0}{\overset{\pi }{{\displaystyle \int }}}(\frac{d^{2}v}{d{y}^2}\left(\frac{1+\cos y}{2}\right))\cos \left(y+\frac{1}{2}\right)y\cos \left(\frac{y}{2}\right)dy$$

(38)

Integrating Equation (38), we obtain

$$c_n=\frac{1}{4\pi }\underset{0}{\overset{\pi }{{\displaystyle \int }}}(\frac{d^{4}v}{d{y}^4}\left(\frac{1+\cos y}{2}\right)){\xi }_n\left(y\right)dy$$

(39)

where ${\xi }_n\left(y\right)=\sin y\left[\frac{1}{n}\left(\frac{\sin \left(n-1\right)y}{n-1}-\frac{\sin \left(n+1\right)y}{n+1}\right)+\frac{1}{n+1}\left(\frac{\sin ny}{n}-\frac{\sin \left(n+2\right)y}{n+2}\right)\right]$.

We now observe that

$$\begin{gathered} \left|{\xi }_n\left(y\right)\right|=\left|\sin y\left[\frac{1}{n}\left(\frac{\sin \left(n-1\right)y}{n-1}-\frac{\sin \left(n+1\right)y}{n+1}\right)+\frac{1}{n+1}\left(\frac{\sin ny}{n}-\frac{\sin \left(n+2\right)y}{n+2}\right)\right]\right| \\ =\left|\left[\frac{1}{n}\left(\frac{1}{n-1}-\frac{1}{n+1}\right)+\frac{1}{n+1}\left(\frac{1}{n}-\frac{1}{n+2}\right)\right]\right| \\ \le \left|\left[\frac{1}{n}\left(\frac{1}{n-1}+\frac{1}{n+1}\right)+\frac{1}{n+1}\left(\frac{1}{n}+\frac{1}{n+2}\right)\right]\right|<\frac{1}{n^2} \end{gathered}$$

(40)

Furthermore, from Equation (40), we can write

$$\left|c_n\right|=\left|\frac{1}{4\pi }\underset{0}{\overset{\pi }{{\displaystyle \int }}}\left(\frac{d^{4}x}{d{y}^4}\left(\frac{1+\cos y}{2}\right)\right){\xi }_n\left(y\right)dy\right|$$

(41)

Thus, from the given condition, we obtain

$$\left|c_n\right|\le \frac{F}{4\pi }{{\displaystyle \int }}_0^{\pi }\left|{\xi }_n\left(y\right)\right|dy$$

which, in view of Equation (40), yields

$$\left|c_n\right|<\frac{M}{\pi n^2}\underset{0}{\overset{\pi }{{\displaystyle \int }}}dy$$

that is,

$$\left|c_n\right|<\frac{M}{n^2}$$

for a given positive constant M.

So, in Equation (36), if we proceed to the limit as$n\to \infty$, we obtain

$${\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n\to \frac{d^{2}v\left(t\right)}{d{t}^2}$$

□

5. Error Estimation

Theorem 3.

Let$\frac{d^{2}v\left(t\right)}{d{t}^2}\in L_{w\left(t\right)}^2\left[0,1\right]$and${\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n$be the$n^{th}$approximation of it obtained by$\left(n+1\right)$shifted Chebyshev polynomials of the third kind. If$\left|\frac{d^{4}v\left(t\right)}{d{t}^4}\right|<M$for some positive constant M, then the following error estimate holds true:

$$E_{n, \frac{d^{2}v\left(t\right)}{d{t}^2}}^2\le \frac{\pi M^2}{12}{G}_3\left(n+1\right)$$

(42)

where

$$E_{n, \frac{d^{2}v\left(t\right)}{d{t}^2}}={\left(\underset{0}{\overset{1}{{\displaystyle \int }}}{\left|\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)-{\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n\right|}^{2}w\left(t\right)dt\right)}^{\frac{1}{2}}$$

and$G_n\left(y\right)$is a Poly-Gamma function given by

$$G_n\left(y\right)={\left(-1\right)}^{n+1}n!{\displaystyle \sum }_{j=0}^{\infty }\frac{1}{{\left(y+j\right)}^{n+1}}$$

Proof. 

The error term is given by

$$\begin{gathered} E_{n, \frac{d^{2}v\left(t\right)}{d{t}^2}}={\left(\underset{0}{\overset{1}{{\displaystyle \int }}}{\left|\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)-{\left(\frac{d^{2}v\left(t\right)}{d{t}^2}\right)}_n\right|}^{2}w\left(t\right)dt\right)}^{\frac{1}{2}} \\ ={\left(\underset{0}{\overset{1}{{\displaystyle \int }}}{\left|{\displaystyle \sum }_{i=n+1}^{\infty }{c}_i{V}_i\left(t\right)\right|}^{2}w\left(t\right)dt\right)}^{\frac{1}{2}} \end{gathered}$$

(43)

By using the orthogonality property of the third-kind Chebyshev polynomials in Equation (43), we obtain

$$E_{n, \frac{d^{2}v\left(t\right)}{d{t}^2}}={\left(\frac{\pi }{2}{{\displaystyle \sum }}_{i=n+1}^{\infty }{\left|c_i\right|}^2\right)}^{\frac{1}{2}}$$

and

$$E_{n, \frac{d^{2}v\left(t\right)}{d{t}^2}}^2\le \frac{\pi P^2}{2}{{\displaystyle \sum }}_{i=n+1}^{\infty }\frac{1}{i^4},$$

(44)

where

$$\left|\frac{d^{4}v\left(t\right)}{d{t}^4}\right|<M$$

for some positive constant M.

Since

$${{\displaystyle \sum }}_{i=n+1}^{\infty }\frac{1}{i^4}=\frac{1}{6}{G}_3\left(n+1\right)$$

(45)

we find from Equations (32) and (33) that

$$E_{n, \frac{d^{2}v\left(t\right)}{d{t}^2}}^2\le \frac{\pi M^2}{12}{G}_3\left(n+1\right)$$

(46)

□

6. Numerical Simulations and Discussion

In this section, numerical simulations of our results are presented by using illustrative figures and numerical tables. We show the effect of the parameters $\mu , \alpha$ and $\rho$ on the concentration of the substance.

6.1. Results for the Spherical Catalyst Model

We first show the effect of the parameters $\mu , \alpha$ and $\rho$ on the concentration of the substance in a spherical catalyst model. Figure 1 illustrates the behaviour of the concentration v(t; α, μ, ρ) of the substance for different values of $\alpha$ at $\mu =1$ and $\rho =1$.

Figure 2 illustrates the behaviour of the concentration $v\left(t; \alpha ,\mu ,\rho \right)$ of the substance for different values of $\rho$ at $\mu =1$ and $\alpha =1$. From Figure 1 and Figure 2, we can see that the concentration of the substance shows a continuous behaviour with the values of the parameters $\alpha$ and $\rho$.

In

Table 1. Comparison of the results with different values of the parameter $\alpha$.

$\mathit{t}$	Present Method $\left(\mathit{\alpha }=0.5\right)$	Method in [16] $\left(\mathit{\alpha }=0.5\right)$	Present Method $\left(\mathit{\alpha }=1\right)$	Method in [16] $\left(\mathit{\alpha }=1\right)$	Present Method $\left(\mathit{\alpha }=1.5\right)$	Method in [16] $\left(\mathit{\alpha }=1.5\right)$
0.0	0.83946	0.83946	0.82874	0.82874	0.81965	0.81965
0.1	0.84301	0.84100	0.83376	0.83040	0.82499	0.82143
0.2	0.84901	0.84562	0.84102	0.83539	0.83276	0.82676
0.3	0.85710	0.85336	0.84995	0.84372	0.84231	0.83565
0.4	0.86801	0.86425	0.86171	0.85542	0.85487	0.84813
0.5	0.88193	0.87834	0.87656	0.87055	0.87069	0.86421
0.6	0.89893	0.89572	0.89458	0.88916	0.88980	0.88393
0.7	0.91917	0.91646	0.91591	0.91132	0.91234	0.90732
0.8	0.94270	0.94067	0.94058	0.93711	0.93825	0.93443
0.9	0.96960	0.96847	0.96858	0.96663	0.96747	0.96530
1.0	1.00000	1.00000	1.00000	1.00000	1.00000	1.00000

, the results are compared between those derived here by the proposed method and by the method which was used in [16] for different choices of$\alpha$ and $\mu =1$ and $\rho =1$.

In

Table 2. Comparison of the results with different values of the parameter $\rho$.

$\mathit{t}$	Present Method $\left(\mathit{\rho }=0.5\right)$	Method in [16] $\left(\mathit{\rho }=0.5\right)$	Present Method $\left(\mathit{\rho }=1.5\right)$	Method in [16]$\left(\mathit{\rho }=1.5\right)$
0.0	0.95541	0.95541	0.64719	0.64719
0.1	0.95708	0.95586	0.65480	0.65039
0.2	0.95924	0.95719	0.66754	0.66003
0.3	0.96164	0.95940	0.68476	0.67623
0.4	0.96473	0.96251	0.70812	0.69918
0.5	0.96858	0.96651	0.73814	0.72918
0.6	0.97322	0.97141	0.77515	0.76657
0.7	0.97868	0.97720	0.81960	0.81179
0.8	0.98496	0.98389	0.87175	0.86535
0.9	0.99206	0.99149	0.93181	0.92786
1.0	1.00000	1.00000	1.00000	1.00000

below, the results are compared between those derived here by our proposed method and by the method used in [16] for different choices of$\rho$ and $\mu =1$ and $\alpha =1$.

In

Table 3. Comparison of the effectiveness factor $\tau$ with different choices of the parameters.

$\mathit{\beta }$	$\mathit{\alpha }$	$\mathit{\rho }$	$\mathit{\mu }$	Present Method	Method in [16]
0.8394623	0.5	1	1	0.9656645	1.003207
0.8287498	1.0	1 1	1 0	0.9943650	1.059701
0.8196580	1.5	1	1	1.0253210	1.098859
0.9554170	1	0.5	1	1.0041781	1.075815
0.8287498	1	1.0	1	0.9943650	1.059701
0.6471921	1	1.5	1	0.9637729	1.029147

below, the effectiveness factor $\left(\tau \right)$ for different choices of the parameters are listed. From Table 3, we can see that $\tau$ increases with the increase in $\alpha$ and decreases with the increase in $\rho$.

In

Table 4. CPU time taken in computation for different values of $n$ at $\beta =0.9554170509$, $\alpha =1, \mu =1$ and $\rho =0.5$.

$\mathit{n}$	Time (s)
4	8.201
7	18.034
10	60.819

, we have listed CPU time of computations for different dimensions of basis.

6.2. Results for the Spherical Biocatalyst Model

We show here the effect of the parameters $\mu , \alpha$ and $\rho$ on the concentration of the substance in a spherical biocatalyst model. Figure 1 shows the behaviour of the concentration $v\left(t; \alpha ,\mu ,\rho \right)$ of the substance for different values of $\rho$ at $\mu =2$. From Figure 3, we can see that the concentration of the substance shows a continuous behaviour with the values of $\rho$.

In

Table 5. Comparison of the results with different choices of the parameter $\rho$.

$\mathit{t}$	Present Method $\left(\mathit{\rho }=1\right)$	Method in [16] $\left(\mathit{\rho }=1\right)$	Present Method $\left(\mathit{\rho }=2\right)$	Method in [16] $\left(\mathit{\rho }=2\right)$	Present Method $\left(\mathit{\rho }=3\right)$	Method in [16] $\left(\mathit{\rho }=3\right)$
0.0	0.79644	0.79644	0.48206	0.48206	0.21426	0.21426
0.1	0.83673	0.80322	0.50437	0.49965	0.21819	0.23594
0.2	0.86367	0.80322	0.50437	0.49965	0.21819	0.23594
0.3	0.84955	0.81262	0.52743	0.52202	0.24253	0.26447
0.4	0.86303	0.82544	0.55976	0.55404	0.28121	0.30696
0.5	0.87859	0.84227	0.60219	0.59640	0.33626	0.36575
0.6	0.89653	0.86338	0.65556	0.64999	0.41109	0.44360
0.7	0.91761	0.88914	0.72092	0.71586	0.50969	0.54338
0.8	0.94184	0.92005	0.79932	0.79522	0.63704	0.66792
0.9	0.96914	0.95674	0.89191	0.88944	0.79866	0.81958
1.0	1.00000	1.00000	1.00000	1.00000	1.00000	1.00000

below, the results are compared between those derived here by the proposed method and those by the method used in [16] for different choices of $\rho$ and $\mu =2$.

In

Table 6. Comparison of the effectiveness factor $\tau$ with different choices of the parameters.

$\mathit{\beta }$	$\mathit{\rho }$	$\mathit{\mu }$	Present Method	Method in [16]
0.7964472	1	1	0.9848465	1.4057760
0.4820697	2	1	0.9468011	0.8943578
0.2142606	3	1	0.7432604	0.6501860

, the $\tau$ values for different choices of parameters are listed. From Table 6, we can see that $\tau$ decreases with the increase in $\rho$.

7. Conclusions

In this article, we have studied the behaviour of the concentration of the pellets at the spherical origin $t=0$. The main advantages of the technique, which we have proposed and used here, are that it handles the singularity at the origin remarkably easily and it provides an accurate solution at the origin and in the neighbourhood of the origin. Our proposed method is easy to handle because both of the mathematical models, which we have considered in this article, are converted into a system of non-linear algebraic equations. From the illustrative figures and numerical tables, the accuracy and efficiency of the proposed method are sufficiently clear. It is also observed that the desired accuracy is achieved by using a lower number of the basis elements. In our future communications, we can use a different class of polynomials and wavelets in order to possibly achieve a better accuracy.