Generalized Shifted Airfoil Polynomials of the Second Kind to Solve a Class of Singular Electrohydrodynamic Fluid Model of Fractional Order

Abstract

In this manuscript, we find the numerical solutions of a class of fractional-order differential equations with singularity and strong nonlinearity pertaining to electrohydrodynamic flow in a circular cylindrical conduit. The nonlinearity of the underlying model is removed by the quasilinearization method (QLM) and we obtain a family of linearized equations. By making use of the generalized shifted airfoil polynomials of the second kind (SAPSK) together with some appropriate collocation points as the roots of SAPSK, we arrive at an algebraic system of linear equations to be solved in an iterative manner. The error analysis and convergence properties of the SAPSK are established in the $L_2$ and $L_{\infty }$ norms. Through numerical simulations, it is shown that the proposed hybrid QLM-SAPSK approach is not only capable of tackling the inherit singularity at the origin, but also produces effective numerical solutions to the model problem with different nonlinearity parameters and two fractional order derivatives. The accuracy of the present technique is checked via the technique of residual error functions. The QLM-SAPSK technique is simple and efficient for solving the underlying electrohydrodynamic flow model. The computational outcomes are accurate in comparison with those of numerical values reported in the literature.

1. Introduction

The present paper develops a (novel) hybrid and efficient computational procedure to deal with the following class of strongly nonlinear and singular boundary value problems (BVPs) with multi-order fractional derivatives given by [1,2]

$$\left\{\begin{array}{c}{\displaystyle {}^{\mathrm{LC}}{D}_p^{\theta }\Phi \left(p\right)+\frac{1}{p}{}^{\mathrm{LC}}{D}_p^{\lambda }\Phi \left(p\right)+{\mathcal{H}}^2\left(1-\frac{\Phi \left(p\right)}{1-\delta \Phi \left(p\right)}\right)=0,}\hfill \\ {\Phi }^{\prime }\left(0\right)=0,\Phi \left(1\right)=0,\hfill \end{array}\right.p\in (0,1),$$

(1)

where $\theta \in (1,2]$ and $\lambda \in (0,1]$ are the fractional orders of derivative operators ${}^{\mathrm{LC}}{D}_p^{\theta },{}^{\mathrm{LC}}{D}_p^{\lambda }$ are described in the sense of Liouville–Caputo. Here, $\Phi \left(p\right)$ stands for the velocity of fluid; p denotes the radial distance from the center of cylindrical conduit; the parameter ${\mathcal{H}}^2$ is the Hartman electric number; and finally, the degree of the nonlinearity of the (1) is denoted by $\delta$. To be more precise, two parameters ${\mathcal{H}}^2$, $\delta$ are defined by [3]

$$\begin{array}{c}\hfill {\mathcal{H}}^2:=\sqrt{\frac{r^2{j}_0}{\mu E_0{m}^2}},\delta :=\frac{m}{j_0}\nabla ·P-1,\end{array}$$

where r is the radius of the cylindrical conduit, $j_0$ denotes the uniform electrical current density at the inlet, the viscosity of the fluid is represented by $\mu$, m signifies the mobility of ion, and $E_0$ represents the electrical field. Finally, by $\nabla ·P$, we denote the pressure gradient, which is assumed to be constant.

Traditionally, for $\theta =2$ and $\lambda =1$, the model problem (1) has been called the electrohydrodynamic (EHD) flow model. It was used to model the motion of ionized fluid particles and their interrelationships with electric field and the surrounding fluid that controls the transport phenomenon in the fluid flow [4]. The issues of existence and uniqueness of the integer-order model are addressed in [5]. It is further shown that this unique solution is monotonically decreasing on $(0,1)$ and is bounded by $1/(1+\delta )$. The non-fractional model has been considered by many researchers due to its vast applications. Numerous numerical and analytical procedures have been developed for $\theta =2$ and $\lambda =1$. The homotopy analysis and optimal homotopy asymptotic approaches were investigated in [6,7,8]. Additionally, the spectral homotopy analysis method was presented in [9]. The authors in [10] devised the pseudo-spectral collocation strategy for solving the EHD model. The least-square method was proposed in [11]. This problem was solved by utilizing the Runge–Kutta (fourth-order), collocation and Galerkin computational methods in [12]. The author of [13] studied the application of (orthonormal) Bernstein functions for the integer-order EHD flow equation. The discrete Adomian decomposition procedure was developed in [14] for EHD flow model. The higher-order B-spline methods based on uniform and ononuniform meshes were investigated in [15,16]. A semi-analytical technique based on Green’s function along with Picard’s and Mann’s fixed-point iteration were developed recently in [17]. Additionally, the Haar wavelet solutions was obtained in [18]. Two other wavelet solutions based on the Bernoulli and Jacobi polynomials were obtained in [19]. Ultimately, an efficient and highly accurate numerical approach based on the Chebyshev polynomial of the second kind was examined in [20].

However, there are only a few studies that have been devoted to fractional-order EHD flow model (1) with two factional parameters $\theta \in (1,2]$ and $\lambda \in (0,1]$. Nowadays, fractional derivatives are widely utilized as a powerful instrument for handling nonlinear real-world phenomena in various disciplines of science and engineering; see [21,22]. Indeed, they help us to understand the nature of the models more deeply than the integer-order counterparts. Let us, however, emphasize that the exact solutions of (most) fractional differential equations do not exist analytically. Therefore, many researchers attempt to find the approximate solutions of differential equations of fractional order. Let us mention the related available numerical methods for (1). The reproducing kernel scheme was first investigated in [1] and recently reconsidered in [23]. The generalized differential transform procedure designed in [2]. The spectral collocation strategies based on three different bases such as Legendre, Chebyshev, and Jacobi polynomials considered in [24] for computing the solutions of (1). Moreover, the analytical solution that relied on a power series formula was constructed in [25]. Finally, the Galerkin and spectral collocation based on the Lucas functions was employed in [4].

Our focus for this manuscript is solving the fractional-order counterpart of EHD flow model numerically. To this end, a hybrid spectral collocation technique is devised using a (generalized) shifted version of the airfoil polynomials of the second kind (SAPSK). Indeed, we first applied the quasilinearization method (QLM) to the original nonlinear EHD flow equation to arrive at a family of linearlized singular equations. Hence, the spectral SAPSK collocation matrix scheme is utilized to solve the aforementioned family of quasilinear equations. It should be stressed that the spectral collocation procedures have been successfully employed to deal diverse interesting mathematical models computationally by utilizing of different basis functions such as clique, Fibonacci, Morgan-Voyce, Genocci, Bernoulli, Jacobi, Legendre, Pell, and Vieta–Lucas polynomials; see references [26,27,28,29,30,31,32,33,34,35,36] for more details. The proposed hybrid QLM-SAPSK solution technique leads to a (family of linear) system of algebraic equations that determine the expansion functions in polynomial form. This is the main advantage of our QLM-SAPSK technique in comparison with other spectral collocation methods developed previously in [4,24] for (1). The other benefit is the order of convergence, which shows the higher-order accuracy of the QLM-SAPSK. We also establish the error analysis of the utilized airfoil basis functions compared with the Lucas and Chebyshev, and Legendre polynomials in [4,24].

The organization of this manuscript is given as follows. A review of fractional calculus is described in Section 2. Some facts on airfoil polynomials of the second kind is illustrated in Section 3. The convergence analysis of SAPSK in both $L_2$ and $L_{\infty }$ is established rigorously. In Section 4, the details of the hybrid QLM-SAPSK technique based on the quasilinearization and collocation technique is provided. Through the computation of the residual error function, the accuracy of the proposed techniques is further estimated in Section 4. Several numerical simulations using various model parameters are carried out in Section 5. The performance of QLM-SAPSK is validated through comparisons with the outcomes of the available existing computational schemes. The conclusion of the study is summarized in Section 6.

2. Liouville–Caputo Fractional Derivative

We provide the fundamental and basic facts from fractional calculus theory, which is utilized in the following sections. For more information and applications, we refer interested readers to the monograph [22] or a recently published expository article [37].

Let us consider $\Gamma (·)$ as the classical Gamma function. We continue by recalling that the fractional integral operator is of the Riemann-Liouville type. It is of order $\theta >0$, defined by

$$\begin{array}{c}\hfill {}_0{\mathcal{I}}_p^{\theta }s\left(p\right):=\frac{1}{\Gamma \left(\theta \right)}{\int }_0^p\frac{s\left(r\right)}{{(p-r)}^{1-\theta }}dr,\end{array}$$

where we assumed $s\left(p\right)\in C_{\kappa }$, $\kappa >-1$. A given real function $s\left(p\right)$, $p>0$ is in the space $C_{\kappa }$, $\kappa \in \mathbb{R}$ if there exits a real number $\gamma \in \mathbb{R}$ and a function $h\left(p\right)\in C^{\infty }([0,\infty )$ such that $s\left(p\right)=p^{\gamma }h\left(p\right)$. Moreover, for a $n\in \mathbb{N}$, we call that $s\left(p\right)\in C_{\kappa }^n$ if and only if $s^{\left(n\right)}\left(p\right)\in C_{\kappa }$.

The definition of fractional Liouville–Caputo derivative is given next.

Definition 1.

Suppose that, for a given $n\in \mathbb{N}$, we have $s\in C_{-1}^n$ and $n-1<\theta <n$. The fractional derivative of the Liouville–Caputo type for the function $s\left(p\right)$ of order θ is defined by

$${}^{\mathrm{LC}}{D}_p^{\theta }s\left(p\right):={}_0{\mathcal{I}}_p^{n-\theta }{D}^{n}s\left(p\right)=\frac{1}{\Gamma (n-\theta )}{\int }_0^p{(p-r)}^{n-\theta -1}{s}^{\left(n\right)}\left(r\right)dr,p>0,$$

where $D=\frac{d}{dp}.$

It is remarked that the derivative operator ${}^{\mathrm{LC}}{D}_p^{\theta }$ is linear. For a constant number C, one obtains the following:

$${}^{\mathrm{LC}}{D}_p^{\theta }C=0.$$

(2)

Let we have the function $s\left(p\right)=p^x$. With the aid of the following property, one can calculate its Liouville–Caputo fractional derivative as follows:

$${}^{\mathrm{LC}}{D}_p^{\theta }{p}^x=\left\{\begin{array}{cc}0,\hfill & \mathrm{for}x\in {\mathbb{N}}_0\mathrm{and}x<\lceil \theta \rceil ,\hfill \\ {\displaystyle \frac{\Gamma (x+1)}{\Gamma (x+1-\theta )}{p}^{x-\theta },}\hfill & \mathrm{for}x\in {\mathbb{N}}_0\mathrm{and}x\ge \lceil \theta \rceil \mathrm{or}x\notin \mathbb{N}\mathrm{and}x>\lfloor \theta \rfloor .\hfill \end{array}\right.$$

(3)

Here, ${\mathbb{N}}_0:=\mathbb{N}\cup \left\{0\right\}$ and $\lfloor ·\rfloor$ denotes the floor function, which gives us the largest integer number equal or less than $\theta$. Additionally, $\lceil ·\rceil$ represents the ceil function that produces the smallest integer number equal or greater than $\theta$.

3. The Shifted Airfoil Polynomials and Their Convergence Results

We first review the definition of the original airfoil polynomials of the second kind. The shifted version of these polynomial is then introduced. In particular, we investigate the convergence analysis of shifted airfoil polynomials in detail below.

3.1. The Airfoil Polynomials: A Shifted Version

The pressure on an airfoil in steady and unsteady subsonic flow is computed by means of the so-called airfoil polynomials. For an extensive list of formulae related to these polynomials, we refer the readers to the monograph [38]. Alternatively, these polynomials have also been named as the Chebyshev of fourth kind; see [39,40].

The definition of the airfoil polynomials of the second kind is given by

$${\mathcal{A}}_r\left(x\right)=\frac{sin\left[(r+\frac{1}{2})\psi \right]}{sin\frac{\psi }{2}},x=cos\psi ,$$

(4)

for $-1\le x\le 1$. Clearly, we have ${\mathcal{A}}_0\left(x\right)=1$. For $r=1$, one may use the trigonometric relation $sin3\psi =4sin\psi {cos}^2\psi -sin\psi$ to obtain ${\mathcal{A}}_1\left(x\right)=1+2x$. The following recursive identity will be used to generate the remaining airfoil polynomials as

$${\mathcal{A}}_{r+1}\left(x\right)=2x{\mathcal{A}}_r\left(x\right)-{\mathcal{A}}_{r-1}\left(x\right),r=1,2,\dots .$$

(5)

The next two polynomials ${\mathcal{A}}_2\left(x\right)$ and ${\mathcal{A}}_3\left(x\right)$ via (5) are obtained as

$${\mathcal{A}}_2\left(x\right)=-1+2x+4x^2,{\mathcal{A}}_3\left(x\right)=8x^3+4x^2-4x-1.$$

Based on these observations, the special values at $x=0,\pm 1$ are obtained as follows:

$$\begin{array}{c}\hfill {\mathcal{A}}_r(-1)={(-1)}^r,{\mathcal{A}}_r\left(0\right)=\pm 1,{\mathcal{A}}_r\left(1\right)=2r+1.\end{array}$$

In the Sturm–Liouville form, they satisfy the following ordinary differential equation

$$\frac{d}{dx}\left[(1-x^2)\omega \left(x\right){\mathcal{A}}_r^{\prime }\left(x\right)\right]=-r(r+1)\omega \left(x\right){\mathcal{A}}_r\left(x\right),r\in \mathbb{N}.$$

(6)

where $\omega \left(x\right):={\left(\frac{1-x}{1+x}\right)}^{\frac{1}{2}}$ stands as the weight function. By exploiting the former relation (6), it is not a difficult task to prove that the family of airfoil polynomials becomes orthogonal with regard to $\omega \left(x\right)$ on $(-1,1)$. To be precise, we have

$${\int }_{-1}^{+1}{\mathcal{A}}_r\left(x\right){\mathcal{A}}_s\left(x\right)\omega \left(x\right)dx=\left\{\begin{array}{cc}\pi ,\hfill & r=s,\hfill \\ 0,\hfill & r\ne s.\hfill \end{array}\right.$$

(7)

The explicit expansion form of the airfoil polynomials of the second kind (APSK) is given by

$${\mathcal{A}}_r\left(x\right)=\frac{1}{2^r}\sum _{\ell =0}^r{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{2r+1}{2\ell +1}\right){(1-x)}^{\ell }{(1+x)}^{r-\ell },r\in {\mathbb{N}}_0.$$

(8)

One can easily obtain the zeros of APSK by solving the equation ${\mathcal{A}}_r\left(x_j\right)=0$. The roots of ${\mathcal{A}}_r\left(x\right)$ of degree r are all distinct and real on $(-1,1)$, given by

$$\begin{array}{c}\hfill x_j=-cos\left[\frac{(2j-1)\pi }{2r+1}\right],j=1,2,\dots ,r.\end{array}$$

(9)

In this work, our goal is to employ the APSK on the unit interval $[0,1]$. Thus, the shifted airfoil polynomials are defined by the change in variable $x=2z-1$, as follows.

Definition 2.

Let ${\mathcal{A}}_r^{\star }\left(z\right)$ represent the shifted airfoil polynomials on $[0,1]$ of degree r. It is defined via the following relation:

$${\mathcal{A}}_r^{\star }\left(z\right)={\mathcal{A}}_r(2z-1),z\in [0,1].$$

(10)

By using the aforesaid transformation in (8), we derive the following explicit form:

$${\mathcal{A}}_r^{\star }\left(z\right)=\sum _{\ell =0}^r{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{2r+1}{2\ell +1}\right){(1-z)}^{\ell }{z}^{r-\ell },r\in {\mathbb{N}}_0.$$

(11)

In accordance with (7), one can obtain a similar orthogonality relation for the shifted APSK ${\left\{{\mathcal{A}}_r^{\star }\left(z\right)\right\}}_{r=0}^{\infty }$. It is straightforward to show that the related weight function is ${\omega }_{\star }\left(z\right)=\sqrt{\frac{1-z}{z}}$ for $z\in (0,1)$. Thus, we obtain

$${\int }_0^1{\mathcal{A}}_r^{\star }\left(z\right){\mathcal{A}}_s^{\star }\left(z\right){\omega }_{\star }\left(z\right)dz=\frac{\pi }{2}{\delta }_{rs},$$

(12)

where ${\delta }_{rs}$ is the well-known Kronecker delta function. We next specify the locations of roots related to shifted APSK. The proof of the next result is straightforward using the relation (9) and Definition 2.

Lemma 1.

The zeros of the shifted airfoil polynomials ${\mathcal{A}}_r^{\star }\left(z\right)$ are within (0,1) and given by

$$z_j=\frac{1}{2}\left(1+x_j\right),j=1,2,\dots ,r,$$

(13)

where $x_j$ are given in (9). These points will be used as the collocation points in the sequel.

Finally, we obtain a modified version of explicit form (11) as a powers of z. To this end, we first employ the following binomial expansion:

$$\begin{array}{c}\hfill {(1-z)}^{\ell }=\sum _{i=0}^{\ell }{(-1)}^{\ell -i}\left(\genfrac{}{}{0pt}{}{\ell }{i}\right)z^{\ell -i}.\end{array}$$

After substituting the forgoing relation into (11) followed by some manipulations, we obtain the following representation form for the shifted airfoil polynomials:

$${\mathcal{A}}_r^{\star }\left(z\right)=\sum _{\ell =0}^{r}c(r,\ell )z^{\ell },r\in {\mathbb{N}}_0,$$

(14)

where

$$\begin{array}{c}\hfill c(r,\ell ):=\sum _{i=r-\ell }^{\ell }{(-1)}^{2i-r+\ell }\left(\genfrac{}{}{0pt}{}{2r+1}{2i+1}\right)\left(\genfrac{}{}{0pt}{}{i}{r-\ell }\right).\end{array}$$

3.2. Convergent and Error Analysis

One of the primary aims of this work is to analyze the convergence properties of the shifted APSK rigorously. First, note that a function $\kappa \left(z\right)\in L_{2,\star }\left([0,1]\right)$ can be stated as a linear of combination of SAPSK. Thus, we can write

$$\kappa \left(z\right)=\sum _{r=0}^{\infty }{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right),z\in [0,1].$$

(15)

In accordance with the orthogonality relation (12), the unknown coefficients ${\delta }_r$ are obtained in the closed-form as

$${\delta }_r:=\frac{2}{\pi }{\int }_0^1{\mathcal{A}}_r^{\star }\left(z\right)\kappa \left(z\right){\omega }_{\star }\left(z\right)dz,r=1,2,\dots .$$

(16)

The following result is used to prove that the expansion series solution (15) is uniformly convergent. For this purpose, we first estimate the coefficients ${\delta }_r$ in (15). Thus, we have

Theorem 1.

Let a function $\kappa \in L_{2,\star }\left([0,1]\right)\cap C^{\left(2\right)}\left([0,1]\right)$ be expressed as (15) and ${\mathcal{M}}_2:={max}_{z\in [0,1]}\left|{\kappa }^{″}\left(z\right)\right|$. Then, the following upper bound for ${\delta }_r$ in (16) is obtained

$$\begin{array}{c}\hfill |{\delta }_r|<c_1{r}^{-4},r>1,\end{array}$$

(17)

where $c_1=\frac{4}{\pi }{\mathcal{M}}_2$.

Proof. 

We proceed by using the substitution $z=\frac{1}{2}+\frac{1}{2}cos\psi =:h\left(\psi \right)$ in (16) to arrive at

$${\delta }_r=\frac{2}{\pi }{\int }_0^{\pi }\kappa \left(h\left(\psi \right)\right)sin[(r+\frac{1}{2})\psi ]sin\frac{\psi }{2}d\psi =\frac{1}{\pi }{\int }_0^{\pi }\kappa \left(h\left(\psi \right)\right)\left\{cos\left[r\psi \right]-cos\left[\right(r+1\left)\psi \right]\right\}d\psi .$$

(18)

By employing the integration by parts on (18) twice, we obtain

$${\delta }_r=\frac{1}{8\pi }{\int }_0^{\pi }{\kappa }^{″}\left(h\left(\psi \right)\right){\omega }_r\left(\psi \right)sin\left(\psi \right)d\psi ,$$

(19)

where

$${\omega }_r\left(\psi \right):=\frac{1}{r}\left(\frac{sin\left(\right(r-1\left)\psi \right)}{r-1}-\frac{sin\left(\right(r+1\left)\psi \right)}{r+1}\right)-\frac{1}{r+1}\left(\frac{sin\left(r\psi \right)}{r}-\frac{sin\left(\right(r+2\left)\psi \right)}{r+2}\right).$$

By our assumption, we know that the second derivative is bounded by ${\mathcal{M}}_2$. Additionally, we always have that $|sin(\psi \left)\right|\le 1$. Thus, we conclude that

$$|{\delta }_r|\le \frac{{\mathcal{M}}_2}{8\pi }\left|{\int }_0^{\pi }{\omega }_r\left(\psi \right)d\psi \right|.$$

(20)

Our task is now to estimate the integral term in (20). We apply the change of variables $p=s\psi$, for $s=r,r\pm 1,r+2$ to obtain

$${\int }_0^{\pi }{\omega }_r\left(\psi \right)d\psi =\frac{1-{(-1)}^{r-1}}{{(r-1)}^{2}r}+\frac{{(-1)}^{r+1}-1}{r{(r+1)}^2}+\frac{{(-1)}^r-1}{r^2(r+1)}+\frac{1-{(-1)}^{r+2}}{(r+1){(r+2)}^2}.$$

By an easy calculation, one shows that for each value of $r>1$, two terms on the right-hand side are zero. Using these facts, we obtain for $r>1$

$${\int }_0^{\pi }{\omega }_r\left(\psi \right)d\psi =\left\{\begin{array}{cc}\frac{2}{(r+1){(r+2)}^2}+\frac{-2}{r^2(r+1)}=\frac{-8}{r^2{(r+2)}^2},\hfill & \mathrm{i}\mathrm{f} r \mathrm{o}\mathrm{d}\mathrm{d},\hfill \\ \frac{2}{{(r-1)}^{2}r}+\frac{-2}{r{(r+1)}^2}=\frac{8}{{(r-1)}^2{(r+1)}^2},\hfill & \mathrm{i}\mathrm{f} r \mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}.\hfill \end{array}\right.$$

With the aid of simple inequalities $r+1<r+2$ and $r-1<r$, we immediately conclude that

$$\left|{\int }_0^{\pi }{\omega }_r\left(\psi \right)d\psi \right|\le \frac{8}{{(r-1)}^2{(r+1)}^2}.$$

Now, by using the inequality $r-1\ge \frac{r}{2}$, which is valid for all $r\ge 2$, we obtain

$$\left|{\int }_0^{\pi }{\omega }_r\left(\psi \right)d\psi \right|<\frac{32}{r^4}.$$

(21)

By placing (21) into (20), we complete the desired result (17). □

We practically consider a cut series solution to approximate $\kappa \left(z\right)$ instead of the infinite series solution given in (15). If (15) is truncated up to its first $(R+1)$ term, then it may be written as

$$\kappa \left(z\right)\approx {\kappa }_R\left(z\right)=\sum _{r=0}^R{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right).$$

(22)

Let ${\kappa }_R$ and ${\kappa }_{R+1}$ be two consecutive approximations of $\kappa \left(z\right)$. By $e_R$, we denote the difference between these two expansion series, given by

$$e_R\left(z\right):={\kappa }_{R+1}\left(z\right)-{\kappa }_R\left(z\right).$$

(23)

Moreover, by ${\parallel g\parallel }_{2,\star }$ we present the weighted $L_{2,\star }$ norm on $[0,1]$ with regard to weight function ${\omega }_{\star }\left(z\right)$. In the weighted $L_2$ norm, an upper bound for the error $e_R$ is obtained in next Theorem.

Theorem 2.

Let suppose that the assumptions of Theorem 1 are valid. Then, the following error estimate holds:

$${\parallel e_R\left(z\right)\parallel }_{2,\star }<c_2{R}^{-4},c_2:=\sqrt{\frac{\pi }{2}}{c}_1=\sqrt{\frac{8}{\pi }}{\mathcal{M}}_2.$$

Here, the constants $c_1,{\mathcal{M}}_2$ are given in (17).

Proof. 

On account of definition of error (23) and using (22), we have

$$\begin{array}{cc}\hfill {\parallel e_R\left(z\right)\parallel }_{2,\star }& ={\parallel {\kappa }_{R+1}\left(z\right)-{\kappa }_R\left(z\right)\parallel }_{2,\star }\hfill \\ \hfill & =\parallel \sum _{r=0}^{R+1}{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right)-\sum _{r=0}^R{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right){\parallel }_{2,\star }\hfill \\ \hfill & ={\parallel {\delta }_{R+1}{\mathcal{A}}_{R+1}^{\star }\left(z\right)\parallel }_{2,\star }\hfill \\ \hfill & =|{\delta }_{R+1}|{\parallel {\mathcal{A}}_{R+1}^{\star }\left(z\right)\parallel }_{2,\star }\hfill \end{array}$$

With the help of orthogonality condition (12) and the result of Theorem 1, we have

$$\begin{array}{cc}\hfill {\parallel e_R\left(z\right)\parallel }_{2,\star }& =|{\delta }_{R+1}|\sqrt{\frac{\pi }{2}}\hfill \\ \hfill & <c_1{(R+1)}^{-4}\sqrt{\frac{\pi }{2}}\hfill \\ \hfill & <c_1\sqrt{\frac{\pi }{2}}{R}^{-4}.\hfill \end{array}$$

□

Our next goal is to find an upper bound for (global) error between the infinite series form in (15) for $\kappa \left(z\right)$ and its truncated series solution ${\kappa }_R\left(z\right)$ in (22). Therefore, we define

$$E_R\left(z\right)=\kappa \left(z\right)-{\kappa }_R\left(z\right).$$

(24)

The first result is devoted to the weighted $L_{2,\star }\left([0,1]\right)$ norm given next.

Theorem 3.

Let suppose that the hypotheses of Theorem 1 hold. Then, the (global) error $E_R\left(z\right)$ in the $L_{2,\star }\left([0,1]\right)$ norm satisfies

$$\begin{array}{c}\hfill {\parallel E_R\parallel }_{2,\star }<c_3{R}^{-\frac{7}{2}},c_3:=\sqrt{\frac{\pi }{14}}{c}_1=\sqrt{\frac{8}{7\pi }}{\mathcal{M}}_2,\end{array}$$

where the constant ${\mathcal{M}}_2$ is given in (17).

Proof. 

By following (15) and (22), we obtain

$$\begin{array}{c}\hfill {\parallel E_R\parallel }_{2,\star }^2=\parallel \sum _{r=0}^{\infty }{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right)-\sum _{r=0}^R{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right){\parallel }_{2,\star }^2=\parallel \sum _{r=R+1}^{\infty }{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right){\parallel }_{2,\star }^2.\end{array}$$

Now, we employ the orthogonality condition (12) to arrive at

$$\begin{array}{c}\hfill {\parallel E_R\parallel }_{2,\star }^2=\frac{\pi }{2}\sum _{r=R+1}^{\infty }{\delta }_r^2.\end{array}$$

It is now sufficient to employ the inequality (17), which was already proved in Theorem 3. Applying it to the preceding relation, one obtains the following:

$$\begin{array}{c}\hfill {\parallel E_R\parallel }_{2,\star }^2\le \frac{\pi }{2}{c}_1^2\sum _{r=R+1}^{\infty }\frac{1}{r^8}.\end{array}$$

(25)

The direct application of the Integral Test from calculus gives us [41]

$$\begin{array}{c}\hfill \sum _{r=R+1}^{\infty }\frac{1}{r^8}\le {\int }_R^{\infty }\frac{ds}{s^8}=\frac{1}{7R^7}.\end{array}$$

To obtain the desired result, we put the forgoing inequality into (25). The proof is finished by taking the square root. □

Still, we are interested in obtaining a bound for the global error (24) in the $L_{\infty }$ norm. To pave the way, the next Lemma will be is first proved:

Lemma 2.

The shifted APSK satisfies the following inequality (for all $r\ge 0$)

$$|{\mathcal{A}}_r^{\star }\left(z\right)|\le 2r+1,\forall z\in [0,1].$$

(26)

Proof. 

In order to prove the result, we first recall that airfoil polynomials of the second kind (4) satisfies [38]

$${\mathcal{A}}_r\left(x\right)=T_r\left(x\right)+(1+x)U_{r-1}\left(x\right),x\in [-1,1].$$

where $T_r\left(x\right)$ and $U_r\left(x\right)$ represent the Chebyshev polynomials of the first and second kinds, respectively. For these classical polynomials, the following relations hold [42]:

$$|T_r\left(x\right)|\le 1,\mathrm{and}|U_r\left(x\right)|\le r+1,\forall \left|x\right|\le 1.$$

By changing variable $x=2z-1$ and then taking the absolute values, we conclude that

$$\begin{array}{c}\hfill |{\mathcal{A}}_r^{\star }\left(z\right)|\le |T_r(2z-1)|+|2z{U}_{r-1}(2z-1)|\le 1+2\left(1\right)\left(r\right)=1+2r,\forall z\in [0,1],\end{array}$$

which completes the proof. □

Theorem 4.

Let suppose that the assumptions of Theorem 1 are valid. Then, the (global) error $E_R\left(z\right)={\sum }_{r=R+1}^{\infty }{\delta }_r{\mathcal{A}}_r^{\star }\left(z\right)$ in the $L_{\infty }\left([0,1]\right)$ norm satisfies

$$\begin{array}{c}\hfill {\parallel E_R\parallel }_{\infty }<c_4{R}^{-2},c_4:=\frac{3}{2}{c}_1=\frac{6}{\pi }{\mathcal{M}}_2,\end{array}$$

(27)

where the constant ${\mathcal{M}}_2$ is given in (17).

Proof. 

By employing the inequality (26) obtained in Lemma (2), we have

$$\begin{array}{cc}\hfill |E_R\left(z\right)|& \le \sum _{r=R+1}^{\infty }|{\delta }_r||{\mathcal{A}}_r^{\star }\left(z\right)|\hfill \\ \hfill & \le \sum _{r=R+1}^{\infty }(2r+1)|{\delta }_r|\le \sum _{r=R+1}^{\infty }3r|{\delta }_r|.\hfill \end{array}$$

We now use the inequality (17) in Theorem 3 to obtain

$$\begin{array}{cc}\hfill |E_R\left(z\right)|& <3c_1\sum _{r=R+1}^{\infty }\frac{1}{r^3}.\hfill \end{array}$$

Again, we utilize the Integral Test [41] for the last inequality to show that

$$\begin{array}{c}\hfill \sum _{r=R+1}^{\infty }\frac{1}{r^3}\le {\int }_R^{\infty }\frac{ds}{s^3}=\frac{1}{2R^2}.\end{array}$$

To arrive at the desired inequality, the supermum is taken over elements of $z\in [0,1]$. □

4. The Hybrid QLM-SAPSK Procedure

Basically, one may apply the spectral collocation scheme based on the SAPSK to the original nonlinear model (1) in a direct manner. However, the main shortcoming of this approach is that we have to solve a non-linear system of equations and it is very time-consuming as the number of SAPSK basis functions, R, is increased. In this respect, the procedure of quasilinearization is adopted to transform the nonlinear EHD flow model (1) into a family of linearized equations. Starting from a rough first approximation, the method converges with quadratic order to the solution of the original problem (1).

In what follows, we describe the fundamental fact about the quasilinearization method (QLM). Once the nonlinear EHD flow model problem (1) is transformed into a family of linear problems, we employ the direct SMV matrix collocation procedure to each linearized equation. Below, this combined technique is called QLM-SAPSK. For more detailed information and a recent applications of QLM, we refer to [43,44,45].

Let us begin by reformulating the nonlinear EHD flow model (1) as

$$L^{\theta ,\lambda }[\Phi ]\left(p\right)=F(p,\Phi \left(p\right)),$$

(28)

where the linear operator $L^{\theta ,\lambda }$ and the nonlinear function F are defined as

$${\displaystyle L^{\theta ,\lambda }[\Phi ]\left(p\right):={}^{\mathrm{LC}}{D}_p^{\theta }\Phi \left(p\right)+\frac{1}{p}{}^{\mathrm{LC}}{D}_p^{\lambda }\Phi \left(p\right),F(p,\Phi \left(p\right)):=-{\mathcal{H}}^2\left(1-\frac{\Phi \left(p\right)}{1-\delta \Phi \left(p\right)}\right).}$$

Assume that ${\Phi }_0\left(p\right)$ is the rough first approximation to $\Phi \left(p\right)$. Then, the process of QLM for (28) can be stated as

$$L^{\theta ,\lambda }\left[{\Phi }_{n+1}\right]\left(p\right)\approx F(p,{\Phi }_n\left(p\right))+F_{\Phi }(p,{\Phi }_n\left(p\right))\left({\Phi }_{n+1}\left(p\right)-{\Phi }_n\left(p\right)\right),n=0,1,\dots .$$

Note that, along with former equations, we have the same initial conditions as in (1).

$$\begin{array}{c}\hfill {\Phi }_{n+1}^{\prime }\left(0\right)=0,{\Phi }_{n+1}\left(1\right)=0.\end{array}$$

(29)

Performing some calculations, the applied QLM gives rise to the following representation for the model (28) as

$$L^{\theta ,\lambda }\left[{\Phi }_{n+1}\right]\left(p\right)-\frac{{\mathcal{H}}^2}{{\left(1-\delta {\Phi }_n\left(p\right)\right)}^2}{\Phi }_{n+1}\left(p\right)=-{\mathcal{H}}^2\left(1+\frac{\delta {\Phi }_n^2\left(z\right)}{{\left(1-\delta {\Phi }_n\left(p\right)\right)}^2}\right),n=0,1,\dots .$$

(30)

To remove the denominator, let us multiply both sides of (30) by ${(1-\delta {\Phi }_n\left(p\right))}^2$ to arrive at

$$\begin{array}{c}\hfill {\left(1-\delta {\Phi }_n\left(p\right)\right)}^2{L}^{\theta ,\lambda }\left[{\Phi }_{n+1}\right]\left(p\right)-{\mathcal{H}}^2{\Phi }_{n+1}\left(p\right)=-{\mathcal{H}}^2\left(1+\delta (1+\delta ){\Phi }_n^2\left(p\right)-2\delta {\Phi }_n\left(p\right)\right).\end{array}$$

(31)

Finally, we introduce the following notations

$$\begin{array}{cc}\hfill {\kappa }_{3,n}\left(p\right)& :={\left(1-\delta {\Phi }_n\left(p\right)\right)}^2,{\kappa }_{2,n}\left(p\right):=\frac{{\kappa }_{3,n}\left(p\right)}{p},{\kappa }_{1,n}\left(p\right):=-{\mathcal{H}}^2,\hfill \\ \hfill {\nu }_n\left(p\right)& :=-{\mathcal{H}}^2\left(1+\delta (1+\delta ){\Phi }_n^2\left(p\right)-2\delta {\Phi }_n\left(p\right)\right),n=0,1,\dots .\hfill \end{array}$$

to write the former Equation (31) in a convenient form as

$$\begin{array}{c}\hfill {\kappa }_{3,n}\left(p\right){}^{\mathrm{LC}}{D}_p^{\theta }{\Phi }_{n+1}\left(p\right)+{\kappa }_{2,n}\left(p\right){}^{\mathrm{LC}}{D}_p^{\lambda }{\Phi }_{n+1}\left(p\right)+{\kappa }_{1,n}\left(p\right){\Phi }_{n+1}\left(p\right)={\nu }_n\left(p\right),n=0,1,\dots .\end{array}$$

(32)

We now are in a position to find the solution of the linearized equations (32) through a matrix collocation procedure relied on the SAPSK. In view of (22), we let the approximate solution of (32) be expressed as a cut series with $(R+1)$ bases as

$${\Phi }_{n+1}\left(p\right)\approx {\chi }_R^{(n+1)}\left(p\right)=\sum _{r=0}^R{\delta }_r^{\left(n\right)}{\mathcal{A}}_r^{\star }\left(p\right),$$

(33)

for $n=0,1,\dots$. We then put the unknown coefficients ${\delta }_r^{\left(n\right)}$ in a vector form as

$${\mathbf{\Delta }}_R^{\left(n\right)}={\left[{\delta }_0^{\left(n\right)}{\delta }_1^{\left(n\right)}\dots {\delta }_R^{\left(n\right)}\right]}^T.$$

By forming the vector of SAPSK

$${\mathit{A}}_R={\left[{\mathcal{A}}_0^{\star }\left(p\right){\mathcal{A}}_1^{\star }\left(p\right)\dots {\mathcal{A}}_R^{\star }\left(p\right)\right]}^T,$$

one can rewrite the approximate solution ${\chi }_R^{(n+1)}\left(p\right)$ in (33) as compact form as

$${\chi }_R^{(n+1)}\left(p\right)={\mathit{A}}_R{\mathbf{\Delta }}_R^{\left(n\right)}.$$

(34)

The following result represents ${\mathit{A}}_R$ as a product of a vector by a constant matrix. The proof of this Lemma is straightforward in view of relation (11) in Definition (2).

Lemma 3.

The vector of SAPSK can be rewritten as

$${\mathit{A}}_R\left(p\right)={\mathit{P}}_R\left(p\right){\mathit{M}}_R,$$

(35)

where ${\mathit{P}}_R\left(p\right)=\left[1p{p}^2\dots p^R\right]$ contains the powers of p, the structured matrix ${\mathit{M}}_R={\left(m_{ij}\right)}_{i,j=0}^R$ is defined by (14), and its elements are given as

$$m_{ij}:=\left\{\begin{array}{cc}0,\hfill & i<j,\hfill \\ c(i,j),\hfill & i\ge j.\hfill \end{array}\right.$$

One can obtain the following corollary by combining two relations (34) and (35) immediately:

Corollary 1.

The approximate solution ${\chi }_R^{(n+1)}\left(p\right)$ in (33) is reformulated as

$${\chi }_R^{(n+1)}\left(p\right)={\mathit{P}}_R\left(p\right){\mathit{M}}_R{\mathbf{\Delta }}_R^{\left(n\right)}.$$

(36)

Remark 1.

When we have the fractional-order derivatives $\theta \in (1,2]$ and $\lambda \in (0,1]$ in (1), it is beneficial to use the generalized SAPSK (GSAPSK) by introducing the local parameter $\alpha >0$. Thus, we obtain ${\mathcal{A}}_r^{\alpha }\left(z\right):={\mathcal{A}}_r^{\star }\left(z^{\alpha }\right)$. This methodology have been successfully applied to some other fractional model problems such as the Brusselator chemical model [26], and the Bratu and Lane–Emden equations [46] by the authors. To be more precise, in our proposed algorithm, we need the monomial vector ${\mathit{P}}_R\left(p\right)$ to be replaced by its generalization:

$${\mathit{P}}_{R,\alpha }\left(p\right):=\left[1p^{\alpha }{p}^{2\alpha }\dots p^{R\alpha }\right].$$

It should be emphasized that the matrix representation (36) is very useful especially if one has to compute the higher-order derivatives of the unknown solutions. For this purpose, we require a calculation the derivatives of the vector ${\mathit{P}}_R\left(p\right)$. In this work, we need to compute the first derivative of the solution. Thus, we have the following:

Proposition 1.

The integer-order derivative of ${\chi }_R^{(n+1)}\left(p\right)$ in (33) is computed as

$$\frac{d}{dp}{\chi }_R^{(n+1)}\left(p\right)={\mathit{P}}_R\left(p\right){\mathit{D}}_R{\mathit{M}}_R{\mathbf{\Delta }}_R^{\left(n\right)},$$

(37)

where the so-called differentiation matrix ${\mathit{D}}_R$ is a zero matrix except for entries $d_{ii+1}=i$ for $i=0,1,\dots ,R$.

Proof. 

We first differentiate (36) with regard to p to obtain

$$\begin{array}{c}\hfill \frac{d}{dp}{\chi }_R^{(n+1)}\left(p\right)=\left(\frac{d}{dp}{\mathit{P}}_R\left(p\right)\right){\mathit{M}}_R{\mathbf{\Delta }}_R^{\left(n\right)}.\end{array}$$

Therefore, our task will be accomplished by just differentiating ${\mathit{P}}_R\left(p\right)$ once. An easy calculation reveals that

$$\begin{array}{c}\hfill \frac{d}{dp}{\mathit{P}}_R\left(p\right)={\mathit{P}}_R\left(p\right){\mathit{D}}_R.\end{array}$$

After combining two former relations, we obtain the desired result. □

In addition to the integer-order derivative of ${\chi }_R^{(n+1)}\left(p\right)$ in (34), one needs the fractional derivatives of order $\theta$ and $\lambda$ as they are in the original Equation (1). To this end, the properties (2) and (3) are essential. Inspired by previous work [46] (Algorithm 1), a simple algorithm can be utilized to compute the $\theta$- and $\lambda$-derivatives of ${\mathit{P}}_R\left(p\right)$. The cost of this algorithm is linear with respect to R. From now on, we let vector forms of these derivatives be

$$\begin{array}{c}\hfill {\mathit{P}}_R^{\theta }\left(p\right):={}^{\mathrm{LC}}{D}_p^{\theta }{\mathit{P}}_R\left(p\right),{\mathit{P}}_R^{\lambda }\left(p\right):={}^{\mathrm{LC}}{D}_p^{\lambda }{\mathit{P}}_R\left(p\right).\end{array}$$

From these notations, we obtain the following results:

Proposition 2.

The fractional-order derivatives of ${}^{\mathrm{LC}}{D}_p^{\gamma }{\chi }_R^{(n+1)}\left(p\right)$ in (33) for $\gamma =\theta ,\lambda$ are calculated as

$${}^{\mathrm{LC}}{D}_p^{\gamma }{\chi }_R^{(n+1)}\left(p\right)={\mathit{P}}_R^{\gamma }\left(p\right){\mathit{M}}_R{\mathbf{\Delta }}_R^{\left(n\right)}.$$

(38)

We are now ready to substitute the SAPSK nodes (13) into the vector representation forms of the unknowns and its derivatives. To do so, we use the roots of ${\mathcal{A}}_{R+1}^{\star }\left(p\right)$ labeled as $p_0,p_1,\dots ,p_R$. By utilizing the following matrix forms:

$$\begin{array}{c}\hfill \mathit{P}:=\left(\begin{array}{c}{\mathit{P}}_R\left(p_0\right)\\ {\mathit{P}}_R\left(p_1\right)\\ ⋮\\ {\mathit{P}}_R\left(p_R\right)\end{array}\right),{\mathit{P}}_{\theta }:=\left(\begin{array}{c}{\mathit{P}}_R^{\theta }\left(p_0\right)\\ {\mathit{P}}_R^{\theta }\left(p_1\right)\\ ⋮\\ {\mathit{P}}_R^{\theta }\left(p_R\right)\end{array}\right),{\mathit{P}}_{\lambda }:=\left(\begin{array}{c}{\mathit{P}}_R^{\lambda }\left(p_0\right)\\ {\mathit{P}}_R^{\lambda }\left(p_1\right)\\ ⋮\\ {\mathit{P}}_R^{\lambda }\left(p_R\right)\end{array}\right),\end{array}$$

(39)

we arrive at the following results:

Lemma 4.

In the matrix formulation, the approximate solution ${\chi }_R^{(n+1)}\left(p\right)$ and its fractional γ-derivatives ${}^{\mathrm{LC}}{D}_p^{\gamma }{\chi }_R^{(n+1)}\left(p\right)$ for $\gamma =\theta ,\lambda$ evaluated at the roots of SAPSK can be rewritten as

$$\begin{array}{cc}\hfill {\mathit{X}}_{n+1}& =\mathit{P}{\mathit{M}}_R{\mathbf{\Delta }}_R^{\left(n\right)},\hfill \\ \hfill {\mathit{X}}_{n+1}^{\gamma }& ={\mathit{P}}_{\gamma }{\mathit{M}}_R{\mathbf{\Delta }}_R^{\left(n\right)},\hfill \end{array}$$

(40)

where the vectors of $\mathit{P}$ and ${\mathit{P}}_{\gamma }$ are already defined by (39) and

$$\begin{array}{c}\hfill {\mathit{X}}_{n+1}:=\left(\begin{array}{c}{\chi }_R^{(n+1)}\left(p_0\right)\\ {\chi }_R^{(n+1)}\left(p_1\right)\\ ⋮\\ {\chi }_R^{(n+1)}\left(p_R\right)\end{array}\right),{\mathit{X}}_{n+1}^{\gamma }:=\left(\begin{array}{c}{}^{\mathrm{LC}}{D}_p^{\gamma }{\chi }_R^{(n+1)}\left(p_0\right)\\ {}^{\mathrm{LC}}{D}_p^{\gamma }{\chi }_R^{(n+1)}\left(p_1\right)\\ ⋮\\ {}^{\mathrm{LC}}{D}_p^{\gamma }{\chi }_R^{(n+1)}\left(p_R\right)\end{array}\right).\end{array}$$

(41)

We now turn to the quasilinear model Equation (32) and place the SAPSK nodes into it. In the matrix formats, we obtain the following representation based on the relations (41) as

$$\begin{array}{c}\hfill {\mathit{K}}_{3,n}{\mathit{X}}_{n+1}^{\theta }+{\mathit{K}}_{2,n}{\mathit{X}}_{n+1}^{\lambda }+{\mathit{K}}_{1,n}{\mathit{X}}_{n+1}={\mathit{V}}_n,\end{array}$$

(42)

where we have used the following notations

$${\mathit{V}}_n:=\left(\begin{array}{c}{\nu }_n\left(p_0\right)\\ {\nu }_n\left(p_1\right)\\ ⋮\\ {\nu }_n\left(p_R\right)\end{array}\right),{\mathit{K}}_{\ell ,n}:=\left(\begin{array}{cccc}{\kappa }_{\ell ,n}\left(p_0\right)& 0& \dots & 0\\ 0& {\kappa }_{\ell ,n}\left(p_1\right)& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\kappa }_{\ell ,n}\left(p_R\right)\end{array}\right),\ell =1,2,3.$$

With the aid of relations (40), one arrives at the following fundamental matrix equation at each iteration n as

$$\begin{array}{c}\hfill \left\{{\mathit{K}}_{3,n}{\mathit{P}}_{\theta }+{\mathit{K}}_{2,n}{\mathit{P}}_{\lambda }+{\mathit{K}}_{1,n}\mathit{P}\right\}{\mathbf{\Delta }}_R^{\left(n\right)}={\mathit{V}}_n.\end{array}$$

(43)

This can be equivalently rephrased as follows:

$${\mathit{Y}}_n{\mathbf{\Delta }}_R^{\left(n\right)}={\mathit{V}}_n,or[{\mathit{Y}}_n;{\mathit{V}}_n],{\mathit{Y}}_n:={\mathit{K}}_{3,n}{\mathit{P}}_{\theta }+{\mathit{K}}_{2,n}{\mathit{P}}_{\lambda }+{\mathit{K}}_{1,n}\mathit{P},n=0,1,\dots .$$

(44)

One still requires the incorporation and implementation of the boundary conditions (29) into the matrix format and then it should be entered into the fundamental matrix Equation (44). We begin by considering the first boundary condition ${\Phi }_{n+1}^{\prime }\left(0\right)=0$. To this end, we use the relation (37) followed by tending p to zero. The resulting matrix form is

$${\mathit{Y}}_n^0{\mathbf{\Delta }}_R^{\left(n\right)}=0,{\mathit{Y}}_n^0:={\mathit{P}}_R\left(0\right){\mathit{D}}_R{\mathit{M}}_R,or[{\mathit{Y}}_n^0;0].$$

For the endpoint boundary condition and in accordance to Corollary 1, it is sufficient to use (36). By approaching $p\to 1$, we obtain the matrix format

$${\mathit{Y}}_n^1{\mathbf{\Delta }}_R^{\left(n\right)}=0,{\mathit{Y}}_n^1:={\mathit{P}}_R\left(1\right){\mathit{M}}_R,or[{\mathit{Y}}_n^1;0].$$

By utilizing these two row matrices $[{\mathit{Y}}_n^0;0]$ and $[{\mathit{Y}}_n^1;0]$, we can replace two rows of the matrix equation $[{\mathit{Y}}_n;{\mathit{V}}_n]$ in (44). Let us used the resulting modified system denoted by

$$[{\widehat{\mathit{Y}}}_n;{\widehat{\mathit{V}}}_n],n=0,1,\dots .$$

(45)

By solving (45), for $r=0,1,\dots ,R$ and $n=1,2,\dots$ in (33), we obtain the unknown coefficients ${\delta }_r^{\left(n\right)}$. It is remarked that, to attain the desired level of accuracy, usually taking $n=5$ is sufficient to solve the system (45) via the QLM-SAPSK technique.

Testing Accuracy via REFs

Typically, closed-form solutions to (1) are practically intractable in particular when we have the fractional orders $1<\theta \le 2$ and $0<\lambda \le 1$. In these cases, usually the technique of residual error (REFs) will help us to compute the accurateness of the present QLM-SAPSK strategy. The key of this technique is that the obtained approximate solution in the series form (33) is substituted into (1) to arrive at the following REFs

$${\displaystyle {\mathcal{R}}_R^{\left(n\right)}\left(p\right):=\left|{}^{\mathrm{LC}}{D}_p^{\theta }{\chi }_R^{\left(n\right)}\left(p\right)+\frac{1}{p}{}^{\mathrm{LC}}{D}_p^{\lambda }{\chi }_R^{\left(n\right)}\left(p\right)+{\mathcal{H}}^2\left(1-\frac{{\chi }_R^{\left(n\right)}\left(p\right)}{1-\delta {\chi }_R^{\left(n\right)}\left(p\right)}\right)\right|\cong 0,}$$

for $n=1,2,\dots$, or we write the REFs more conveniently as

$${\mathcal{R}}_R^{\left(n\right)}\left(p\right):=|\left(1-\delta {\chi }_R^{\left(n\right)}\left(p\right)\right)\left[p{}^{\mathrm{LC}}{D}_p^{\theta }{\chi }_R^{\left(n\right)}\left(p\right)+{}^{\mathrm{LC}}{D}_p^{\lambda }{\chi }_R^{\left(n\right)}\left(p\right)\right]+p{\mathcal{H}}^2\left(1-(1+\delta ){\chi }_R^{\left(n\right)}\left(p\right)\right)|\cong 0,$$

(46)

which is obtained from the former one by multiplying it by $p·\left(1-\delta {\chi }_R^{\left(n\right)}\left(p\right)\right)$. Next, we calculate the maximum values of the residual error norm (for a fixed n) by defining

$${\mathcal{E}}_{\infty }\equiv {\mathcal{E}}_{\infty }^R:=\underset{p\in [0,1]}{max}{\mathcal{R}}_R^{\left(n\right)}\left(p\right).$$

(47)

To show the high-order accuracy of our method and in order to justify our theoretical findings, the numerical order of convergence (${\mathrm{o}\mathrm{r}\mathrm{d}}_R^{\infty }$) is defined via the following relation:

$${\mathrm{o}\mathrm{r}\mathrm{d}}_R^{\infty }:=\frac{ln{\mathcal{E}}_{\infty }^R-ln{\mathcal{E}}_{\infty }^{2R}}{ln2}.$$

(48)

5. Experimental Results and Simulations

The aim would be to apply two proposed matrix algorithm based on the SAPSK to the EHD flow model (1). Numerical simulations and MATLAB plots have been sketched to describe the main outcomes derived in this manuscript. Different values of parameters ($\theta ,\lambda ,{\mathcal{H}}^2,\delta$) are considered to show the performance of the QLM-SAPSK technique. The platform of simulations is Matlab software version 2021a on a personal laptop computer with a 2.2 GHz Intel Core i7-10870H processor and 16 GB of RAM.

The QLM parameter $n=5$ is chosen in the computational results below. In the QLM-SAPSK, we use the initial approximation ${\Phi }_0\left(p\right)=0$ or we take it as ${\Phi }_0\left(p\right)=\frac{1}{1+\delta }$. The latter case is usually employed when we have a large ${\mathcal{H}}^2$. Some comparisons are also performed with the available existing numerical models, i.e., the Lucas spectral collocation method (LSCM) [4], the Galerkin method (GM) based on the Lucas polynomials [4], the least square method (LSM) [11], the cubic B-spline (CBS) function of six order [16], the Haar wavelet (HW) [18] and the discrete Adomian decomposition scheme (DADS) [14], the generalized differential transform method (GDTM) [2], and the reproducing kernel Hilbert space method (RKHSM) [1,2].

Let us consider the integer-order cases. Thus, we set the following parameters in the model EHD flow model (1)

$$\theta =2,\lambda =1,\delta =1,{\mathcal{H}}^2=1.$$

Using the QLM-SAPSK collocation matrix technique with $R=7$, one obtains

$$\begin{array}{cc}\hfill {\chi }_7^{\left(5\right)}\left(p\right)& =-0.0002357p^7+0.00170448p^6-0.000712281p^5-0.0178739p^4-0.000166102p^3\hfill \\ & -0.186131p^2+0.203415.\hfill \end{array}$$

The above approximation is plotted in Figure 1. To show that our proposed approach converges numerically and to justify the theoretical findings, we also use larger values of $R=14,21$. The results of REFs associated to these number of basses are shown in Figure 1, right panel.

We also validate our numerical results obtained by the above parameters by comparing the solutions with the outputs of well-established computational procedures, as shown in

Table 1. A comparison of numerical outcomes/REFs in QLM-SAPSK with $\theta =2,\lambda =1$, $\delta ,{\mathcal{H}}^2=1$, $R=7,14$, and various $p\in [0,1]$.

	QLM-SAPSK				LSM [11]	GM [4]	LSCM [4]
$\mathit{p}$	${\mathbf{\chi }}_{\mathbf{7}}^{\left(\mathbf{5}\right)}\left(\mathit{p}\right)$	${\mathcal{R}}_{\mathbf{7}}^{\left(\mathbf{5}\right)}\left(\mathit{p}\right)$	${\mathbf{\chi }}_{\mathbf{14}}^{\left(\mathbf{5}\right)}\left(\mathit{p}\right)$	${\mathcal{R}}_{\mathbf{14}}^{\left(\mathbf{5}\right)}\left(\mathit{p}\right)$	$\mathit{N}=\mathbf{10}$	${\mathit{N}}_{\mathit{G}}=\mathbf{5}$	${\mathit{N}}_{\mathit{C}}=\mathbf{7}$
$0.0$	$0.20341502$	$0.0000\times {10}^{-0}$	$0.203415795896659$	$0.0000\times {10}^{-00}$	$0.20343243$	$0.20343574$	$0.20342786$
$0.1$	$0.20155175$	$1.8830\times {10}^{-6}$	$0.201552363573421$	$2.3386\times {10}^{-13}$	$0.20156532$	$0.20157001$	$0.20155941$
$0.2$	$0.19593971$	$1.5963\times {10}^{-8}$	$0.195940129562329$	$2.9106\times {10}^{-14}$	$0.19594756$	$0.19594122$	$0.19593886$
$0.3$	$0.18651339$	$3.0982\times {10}^{-7}$	$0.186513678715622$	$5.3721\times {10}^{-14}$	$0.18651760$	$0.18651537$	$0.18654911$
$0.4$	$0.17316508$	$1.5520\times {10}^{-7}$	$0.173165302575071$	$4.7990\times {10}^{-14}$	$0.17316801$	$0.17317053$	$0.17316782$
$0.5$	$0.15574680$	$1.5535\times {10}^{-7}$	$0.155746961727708$	$4.0416\times {10}^{-14}$	$0.15574958$	$0.15574892$	$0.15573556$
$0.6$	$0.13407288$	$1.6136\times {10}^{-7}$	$0.134073001192206$	$3.6251\times {10}^{-14}$	$0.13407547$	$0.13407591$	$0.13408008$
$0.7$	$0.10792349$	$8.2611\times {10}^{-8}$	$0.107923574051674$	$3.9266\times {10}^{-14}$	$0.10792535$	$0.10792460$	$0.10792695$
$0.8$	$0.07704865$	$2.1210\times {10}^{-7}$	$0.077048701459288$	$5.2063\times {10}^{-14}$	$0.07704953$	$0.07704912$	$0.07705274$
$0.9$	$0.04117284$	$1.9348\times {10}^{-7}$	$0.041172865563138$	$3.4394\times {10}^{-14}$	$0.04117309$	$0.04117301$	$0.04118162$

. These approaches are the LSCM [4] with $N_c=7$, the GM based on the Lucas polynomials [4] with $N_G=5$, and the LSM [11] with $N=10$. Our presented results are obtained by using $R=7$ and $R=14$. The resulting REFs are further displayed in this table for justification.

Using the above parameters, we also investigate the behavior of the achieved errors ${\mathcal{E}}_{\infty }$ in the QLM-SAPSK and compute the associated order of convergence, ${\mathrm{o}\mathrm{r}\mathrm{d}}_R^{\infty }$, in

Table 2. The outcomes of ${\mathcal{E}}_{\infty }$ error norms, the related ${\mathrm{o}\mathrm{r}\mathrm{d}}_R^{\infty }$, and the spent CPU time with $\delta ,{\mathcal{H}}^2=1$, $\theta =2,\lambda =1$, and different R.

	QLM-SAPSK				CBS [16]
R	${\mathcal{E}}_{\mathit{\infty }}$	${\mathbf{ord}}_{\mathit{R}}^{\mathit{\infty }}$	CPU(s)	n	MAE	ROC	CPU(s)
2	$5.4729\times {10}^{-1}$	−	$0.49927$	16	$1.7062\times {10}^{-11}$	−	$0.1745$
4	$8.4250\times {10}^{-3}$	$6.0215$	$0.62609$	32	$3.5011\times {10}^{-13}$	$5.6068$	$0.4577$
8	$1.9782\times {10}^{-7}$	$15.378$	$0.89085$	64	$5.8009\times {10}^{-15}$	$5.9154$	$1.0182$
16	$9.6751\times {10}^{-14}$	$20.963$	$1.47224$	128	$9.1434\times {10}^{-17}$	$5.9874$	$2.8513$

. In this respect, we utilize $R=2,4,8,16$, and exploit the relation (48). For each R, the required CPU times (in seconds) are also tabulated in Table 2. It should be remarked that the spent CPU time is for solving the last modified system of equations $[{\widehat{\mathit{Y}}}_n;{\widehat{\mathit{V}}}_n]$ in (45). For comparison, the results of maximum absolute errors (MAE) together with the related rate of convergence (ROC) and the needed CPU time obtained via the cubic B-spline (CBS) function of six order [16] are reported in Table 2. By looking at Table 2, one can clearly observe that our proposed QLM-SAPSK with a lower computational complexity produces more accurate results in comparison with the CBS approach.

In addition to the results presented at Table 2, we also perform a comparison between our outcomes and the results of the spectral collocation procedures based on the classical polynomials such as Jacobi, Legendre, and Chebyshev reported in [24] with $N=10$ bases. The results are shown in

Table 3. A comparison of ${\mathcal{E}}_{\infty }$ error norms utilizing $R=10$, $\delta =0.5$, $\theta =2,\lambda =1$, and different ${\mathcal{H}}^2$.

	QLM-SAPSK	Jacobi [24]	Legendre [24]	Chebyshev [24]	HW [18]	DADS [14]
${\mathcal{H}}^{\mathbf{2}}$	$\mathit{R}=\mathbf{10}$	$\mathit{N}=\mathbf{10}$	$\mathit{N}=\mathbf{10}$	$\mathit{N}=\mathbf{10}$	$\mathit{K}=\mathbf{16}$	$\mathit{N},\mathit{n}=\mathbf{5}$
$0.5$	$2.2246\times {10}^{-11}$	$2.4881\times {10}^{-9}$	$1.3623\times {10}^{-10}$	$6.4622\times {10}^{-10}$	−	$1.6287\times {10}^{-8}$
$1.0$	$9.8648\times {10}^{-10}$	$6.0006\times {10}^{-9}$	$3.1226\times {10}^{-9}$	$3.3842\times {10}^{-9}$	$4.189\times {10}^{-5}$	$9.4029\times {10}^{-7}$
$2.0$	$3.4171\times {10}^{-8}$	$3.3174\times {10}^{-7}$	$1.7631\times {10}^{-7}$	$1.1158\times {10}^{-7}$	$2.421\times {10}^{-7}$	$4.5114\times {10}^{-5}$
$4.0$	$9.9502\times {10}^{-7}$	$8.2384\times {10}^{-6}$	$4.5498\times {10}^{-6}$	$2.9458\times {10}^{-6}$	$6.733\times {10}^{-6}$	$1.4724\times {10}^{-5}$

in which we used a fixed $\delta =0.5$, $R=10$, and ${\mathcal{H}}^2$ varies as $0.5,1,2,4$. Moreover, the outcomes of two other well-established methods, i.e., the Haar wavelet (HW) [18] and the discrete Adomian decomposition scheme (DADS) [14] are further tabulated in Table 3 for comparison. It can be evidently observed that the QLM-SAPSK produces more accurate results than other available approaches.

Figure 2, left plot, shows graphical representations of numerical solutions associated with the results reported in Table 3. Note that the numerical solution with ${\mathcal{H}}^2=10$ is also depicted in this figure. To see the impact of the nonlinearity parameter $\delta$ on the numerical solutions, we also consider $\delta =2$, which is greater than unity. The right picture on Figure 2 presents the obtained results with various parameter ${\mathcal{H}}^2=0.5,1,2,4,10$, but with $\delta =2$. One should emphasize that for plotting the later figure, we have used the initial approximation $1/(1+\delta )$ rather than the zero function. For the larger values of ${\mathcal{H}}^2$, this initial guess is more effective and gives rise to a higher accuracy.

Next, we fix (a large value of) ${\mathcal{H}}^2$ and vary the parameter $\delta =0.1,0.3,0.6,1,2$ to investigate the influence of the nonlinearity on the computations. Firstly, let us take ${\mathcal{H}}^2=2$. Figure 3 shows the results of approximations in accordance with these values of $\delta$. The achieved REFs are further shown on the same figure and on the right part. Obviously, the magnitude of the REFs is an increasing function when we increase the value of $\delta$ as the strength of nonlinearity. In the other words, the quality of the approximation deteriorates as one increases $\delta$. The same situation is also happened if one increases the value of the Hartmann electric number ${\mathcal{H}}^2$. In these cases, the remedy is to increase the number of basis functions accordingly. For ${\mathcal{H}}^2=25$ and ${\mathcal{H}}^2=100$, we plot the approximate solutions related to different values of $\delta$ in Figure 4. Here, we utilized $R=20$ in the computations.

Note that it was shown by Paullet [5] that all approximate solutions are bounded by $\frac{1}{\delta +1}$. This fact can be seen from the plotted Figure 2, Figure 3 and Figure 4. For more discussions on using different large values of the nonlinearity parameters $\delta$ and ${\mathcal{H}}^2$ in the case of integer-order derivatives, we refer to a recent work [20]. Next, we consider the fractional-order cases.

The Fractional-Order Cases

In the second part, our concentration is based on using the fractional-order derivatives $\theta$ and $\lambda$. The aim is to examine the influences of utilizing the simultaneous employment of nonlinearity parameters $\delta ,{\mathcal{H}}^2$ and $\theta ,\lambda$ on the computed approximate solutions. We show these effects through tables and figures. We first pay attention to the following parameters for the multi-order model (1), as previously considered by [1,23]

$$\theta =1.9,\lambda =0.9,\delta =0.5,{\mathcal{H}}^2=1.$$

We first make a comparison between three different choice of the local parameter $\alpha$ as $\lambda ,1,\theta$. By taking $R=10$, the numerical solutions using these parameters and via the QLM-SAPSK/GSAPSK are presented in Figure 5, left picture. Clearly, a good alignment between two approximate solutions obtained by $\alpha =1$ and $\alpha =\theta$ is observed. To choose which $\alpha$ gives us a better resolution, we plot the associated REFs in the same figure and on the right panel. One obviously notices that applying $\alpha =\theta$ yields a smallest absolute values of residual error compared with $\alpha =1$ and $\alpha =\lambda$. On the other hand, by increasing the number of basis functions R, we gain more accurate results.

Therefore, in the next experimental results, we use $R=20$ in the QLM-GSAPSK approach. We fix $\delta =0.5$ and employ various ${\mathcal{H}}^2=1,2,5,10$ in the computations. The results of approximate solutions with $\alpha =\theta$ together with the achieved REFs for the aforesaid parameters are shown in Figure 6. It can be obviously observed that by increasing R, we obtain the desired level of accuracy. In terms of the achieved ${\mathcal{E}}_{\infty }$ error norm, our results presented in Figure 6 can be compared with those graphical plots shown by the generalized differential transform method (GDTM) with $M=140$ terms in [2] (see Figure 3). The maximum absolute errors achieved by the GDTM are $1\times {10}^{-13}$, $4\times {10}^{-11}$, $1.2\times {10}^{-7}$, $5\times {10}^{-5}$ for using ${\mathcal{H}}^2=1,2,5,10$, respectively. It is evident that our results with $R=20$ are more accurate than those obtained via GDTM with $M=140$ terms.

Some precise comparisons are also performed in

Table 4. A comparison of numerical outcomes/REFs in QLM-GSAPSK with $\theta =1.9,\lambda =0.9$, $\delta =0.5,{\mathcal{H}}^2=2$, $R=10$, $\alpha =1,1.9$, and various $p\in [0,1]$.

	QLM-SAPSK
p	$\mathit{\alpha }=1.9$	REFs	$\mathit{\alpha }=1$	REFs	RKHSM [1]	RKHSM [23]
$0.0$	$0.368955980597229$	$0.0000\times {10}^{-00}$	$0.3687331653$	$0.0000\times {10}^{-0}$	$0.381236310$	$0.3825201644$
$0.1$	$0.365172617552478$	$2.6995\times {10}^{-10}$	$0.3650687781$	$5.7028\times {10}^{-5}$	$0.374950080$	$0.3754454339$
$0.2$	$0.354724901744587$	$1.9371\times {10}^{-10}$	$0.3546658303$	$9.9371\times {10}^{-6}$	$0.359968470$	$0.3600500656$
$0.3$	$0.337827276925944$	$1.5404\times {10}^{-11}$	$0.3377883382$	$1.0603\times {10}^{-5}$	$0.342105510$	$0.3434211029$
$0.4$	$0.314292044298967$	$1.2527\times {10}^{-11}$	$0.3142653985$	$9.2000\times {10}^{-6}$	$0.315723980$	$0.3158370163$
$0.5$	$0.283712502798752$	$7.4387\times {10}^{-12}$	$0.2836938281$	$4.7247\times {10}^{-6}$	$0.284128540$	$0.2874089253$
$0.6$	$0.245509789495880$	$4.2012\times {10}^{-12}$	$0.2454968764$	$2.1833\times {10}^{-7}$	$0.244546771$	$0.2475168041$
$0.7$	$0.198955842096497$	$2.7814\times {10}^{-12}$	$0.1989472632$	$3.2885\times {10}^{-6}$	$0.197105564$	$0.2007161568$
$0.8$	$0.143189521840292$	$2.6542\times {10}^{-12}$	$0.1431843235$	$4.3863\times {10}^{-6}$	$0.141053509$	$0.1443651864$
$0.9$	$0.077231456051811$	$1.7987\times {10}^{-12}$	$0.0772290694$	$4.8701\times {10}^{-6}$	$0.075613972$	$0.0777826865$

for the fractional-orders $\theta =1.9$ and $\lambda =0.9$. Here, we utilize $R=10$, $\delta =0.5$, and ${\mathcal{H}}^2=2$. The outcomes of the previously published method, namely the reproducing kernel Hilbert space method (RKHSM) [1,2] are reported in this table for validation. As illustrated before, using the parameter $\alpha$ as the fractional order $\theta$ gives rise to a better accuracy.

It should be emphasized that only one to two digits agreements between our results, and the outcomes of the RKHSM are seen in Table 4. Since the authors in [1,23] did not report the related errors, it is not clear how many digits of accuracy they have. To show that our method delivers more accurate results compared to the previous RKHSM, we consider the case $\alpha =1.9$ and use different $R=5,10,15$. The numerical results are tabulated in

Table 5. A comparison of numerical outcomes/REFs in QLM-GSAPSK with $\theta =1.9,\lambda =0.9$, $\delta =0.5,{\mathcal{H}}^2=2$, $R=5,10,15,20$, $\alpha =1.9$, and various $p\in [0,1]$.

p	$\mathit{R}=5$	$\mathit{R}=10$	$\mathit{R}=15$	$\mathit{R}=20$
$0.0$	$0.368913634324475$	$0.368955980597229$	$0.368955980703425$	$0.368955980703425$
$0.1$	$0.365132468934994$	$0.365172617552478$	$0.365172617645275$	$0.365172617645275$
$0.2$	$0.354690183822644$	$0.354724901744587$	$0.354724901811129$	$0.354724901811129$
$0.3$	$0.337799646966645$	$0.337827276925944$	$0.337827276969954$	$0.337827276969954$
$0.4$	$0.314271598052510$	$0.314292044298967$	$0.314292044328889$	$0.314292044328889$
$0.5$	$0.283698182319713$	$0.283712502798752$	$0.283712502819713$	$0.283712502819713$
$0.6$	$0.245500072203232$	$0.245509789495880$	$0.245509789510364$	$0.245509789510364$
$0.7$	$0.198949413418262$	$0.198955842096497$	$0.198955842106158$	$0.198955842106158$
$0.8$	$0.143185607077043$	$0.143189521840292$	$0.143189521846127$	$0.143189521846127$
$0.9$	$0.077229663096623$	$0.077231456051811$	$0.077231456054508$	$0.077231456054508$

. It should be remarked that the results obtained using $R=15$ and $R=20$ are the same up to 16 digits of accuracy. By looking at Table 5, we infer that by increasing R, the number of repeating decimal digits is increased, which indicates that the proposed QLM-GSAPSK technique is convergent.

Next, we consider the fractional-order parameters as [2]

$$\theta =1.5,\lambda =0.5.$$

By using $R=10$, $\delta =0.5$, and various ${\mathcal{H}}^2=1,2,5,10$, the results of approximations are plotted in Figure 7. As expected, by increasing the nonlinear Hartmann parameter ${\mathcal{H}}^2$, the achieved errors increase. The $L_{\infty }$ error norms (47) along with the associated order of convergences computed via (48) are also tabulated in

Table 6. The outcomes of ${\mathcal{E}}_{\infty }$ error norms and the related ${\mathrm{o}\mathrm{r}\mathrm{d}}_R^{\infty }$ with $\delta =0.5,{\mathcal{H}}^2=1,2,5,10$, $\theta =1.5,\lambda =0.5$, $\alpha =\theta$, and different R.

	${\mathcal{H}}^2=1$		${\mathcal{H}}^2=2$		${\mathcal{H}}^2=5$		${\mathcal{H}}^2=10$
$\mathit{R}$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$
2	$8.2045\times {10}^{-01}$	−	$3.1534\times {10}^{+00}$	−	$1.7334\times {10}^{+01}$	−	$5.6053\times {10}^{+01}$	−
4	$2.6454\times {10}^{-03}$	$8.2768$	$3.4686\times {10}^{-02}$	$6.5064$	$3.0398\times {10}^{-01}$	$5.8335$	$3.6669\times {10}^{+00}$	$3.9345$
8	$2.5453\times {10}^{-08}$	$16.665$	$7.9286\times {10}^{-06}$	$12.095$	$1.5185\times {10}^{-03}$	$7.6452$	$7.4638\times {10}^{-02}$	$5.6181$
16	$8.3563\times {10}^{-16}$	$24.860$	$1.9781\times {10}^{-13}$	$25.256$	$2.7494\times {10}^{-08}$	$15.753$	$5.5263\times {10}^{-06}$	$13.721$

for these values of ${\mathcal{H}}^2$. The high-order accuracy of the proposed QLM-GSAPSK technique applied to the model (1) is evidently seen.

Similarly, we next examine the influence of the nonlinearity term $\delta$ for a fixed ${\mathcal{H}}^2$ in the case of fractional-order derivatives $\theta =1.5,\lambda =0.5$. In this respect, we fix ${\mathcal{H}}^2=0.5$ and vary $\delta =0.5,1,2,5$. The outcomes of error norms ${\mathcal{E}}_{\infty }$ together with related ${\mathrm{o}\mathrm{r}\mathrm{d}}_R^{\infty }$ are displayed in

Table 7. The outcomes of ${\mathcal{E}}_{\infty }$ error norms and the related ${\mathrm{o}\mathrm{r}\mathrm{d}}_R^{\infty }$ with ${\mathcal{H}}^2=0.5,\delta =0.5,1,2,5$, $\theta =1.5,\lambda =0.5$, $\alpha =\theta$, and different R.

	$\mathit{\delta }=0.5$		$\mathit{\delta }=1$		$\mathit{\delta }=2$		$\mathit{\delta }=5$
$\mathit{R}$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$	${\mathcal{E}}_{\infty }$	${\mathbf{ord}}_{\mathit{R}}^{\infty }$
2	$2.0904\times {10}^{-01}$	−	$2.1260\times {10}^{-01}$	−	$2.2009\times {10}^{-01}$	−	$2.4592\times {10}^{-01}$	−
4	$1.7041\times {10}^{-04}$	$10.2606$	$1.3556\times {10}^{-04}$	$10.6150$	$1.0364\times {10}^{-03}$	$7.7304$	$1.8625\times {10}^{-02}$	$3.7229$
8	$1.8749\times {10}^{-11}$	$23.1157$	$3.4860\times {10}^{-09}$	$15.2470$	$1.3152\times {10}^{-07}$	$12.944$	$1.8081\times {10}^{-04}$	$6.6866$
16	$2.7613\times {10}^{-16}$	$16.0511$	$2.7841\times {10}^{-16}$	$23.5779$	$1.7266\times {10}^{-15}$	$26.183$	$7.9406\times {10}^{-09}$	$14.475$

. Clearly, the degree of achieved accuracy decreases if one increases the factor of nonlinearity in the model under consideration.

In the ultimate stage, our aim is to investigate the impact of utilizing diverse values of the factional-order derivatives $\theta$ and $\lambda$ while two other nonlinearity parameters are assumed to be fixed. To this end, we take different pair of orders $(\theta ,\lambda )$ as $(1.2,0.6)$, $(1.4,0.7)$, $(1.6,0.8)$, and $(1.8,0.9)$ and the integer-order $(2,1)$ as a reference. The numerical results using ${\mathcal{H}}^2=0.5$ and a relatively large $\delta =4$ are presented in Figure 8. To use ${\mathcal{H}}^2=50$ and $\delta =0.5$, we utilize $R=30$ in the last experiment. Figure 9 presents the results for different values of $(\theta ,\lambda )$. In both figures, the associated REFs are also depicted.

6. Conclusions

Generalized shifted airfoil polynomials of the second kind (GSAPSK) along with the quasilinearization method (QLM) are utilized to acquire the approximate solutions to a class of fractional-order differential equations with singularity and strong nonlinearity pertaining to electrohydrodynamic (EHD) flow in a circular cylindrical conduit. The fractional operators are considered as the Liouville–Caputo derivative. Besides two fractional orders $\theta$, $\lambda$, the solutions of the EHD flow model depend strongly on the nonlinearity parameters ${\mathcal{H}}^2$ and $\delta$ and, together with the existence of singularities, make handling the underlying model numerically difficult. From a practical point of view, we remove the intrinsic nonlinearity by employing the technique of QLM, which leads to a family of linearized singular equations of fractional order. Then, the GSAPSK collocation matrix technique solved this family using diverse small as well as large model parameters $\theta$, $\lambda$, ${\mathcal{H}}^2$, and $\delta$. In particular, we have shown that the SAPSK expansion series is convergent in the $L_{\infty }$ and weighted $L_2$ norms. Numerical simulations by employing different model parameters are given to show the utility and applicability of the presented QLM-GSAPSK. For validation, our results were compared with those obtained via recently developed analytical and computational techniques in the literature to validate our proposed technique. Overall, the numerical outcomes justify the high accuracy and efficiency of the QLM-GSAPSK method in solving strongly nonlinear EHD flow equations with multi-order derivatives. Furthermore, the achieved numerical orders of convergence indicate that the proposed has an exponential rate of convergence.