A third-order two-step numerical scheme for heat and mass transfer of chemically reactive radiative MHD power-law fluid

Abstract

A two-stage third-order numerical scheme is proposed for solving ordinary differential equations. The scheme is explicit and implicit type in two stages. First, the stability region of the scheme is found when it is applied to the linear equation. Further, the stability conditions of the scheme are found using a linearized homogenous set of differential equations. This set of equations is obtained by applying transformations on the governing equations of heat and mass transfer of incompressible, laminar, steady, two-dimensional, and non-Newtonian power-law fluid flows over a stretching sheet with effects of thermal radiations and chemical reaction. The proposed scheme with an iterative method is employed in two different forms called linearized and non-linearized. But it is found that the non-linearized approach performs better than the linearized one when residuals are compared through plots. Additionally, the proposed scheme is compared to the second-order central finite difference method for second-order non-linear differential equations and the Keller-Box/trapezoidal method for a linear differential equation. It is determined that the proposed scheme is more effective and computationally less expensive than the standard/classical finite difference methods. Moreover, the impact of magnetic parameter, radiation parameter, modified Prandtl and Schmidt numbers for power-law fluid, and chemical reaction rate parameter on velocity, temperature, and concentration profiles are displayed through graphs and discussed. The power-law fluid’s heat and mass transfer simulations are also carried out with varying flow behavior index, sheet velocity, and mass diffusivity. We hoped that this effort would serve as a guide for investigators tasked with resolving unresolved issues in the field of enclosures used in industry and engineering.

Keywords

Third-order scheme stability power-law fluid linearization thermal radiations

Introduction

Numerical methods are considered for solving especially those problems that cannot be solved exactly. The basic type of numerical methods can be classified into two categories. One category is based on using consecutive grid points for ordinary differential equations and consecutive time levels for partial differential equations for time discretizations. These methods can be constructed by using interpolation of polynomials in either an explicit or implicit manner, which are called multistep methods. The other category of numerical methods is based on the use of different stages. The family of Runge-Kutta methods falls into this category. These methods can either be constructed using Butcher’s tableau or the Taylor series are given in this work. In literature, some Runge-Kutta methods have been found which can be used to solve stiff equations. In Kennedy and Carpenter,¹ implicit-explicit fourth and fifth-order Runge-Kutta methods have been given, which used explicit, single diagonally implicit methods, and embedded method was used to control error. Both methods were tested for Van der Pol and Kaps’ problems, and the improvement over existing methods was shown. It was mentioned in Bouhamidi and Jbilou² that implicit Runge-Kutta methods are popular for solving stiff ordinary differential equations. Still, their main drawback is the computation cost for integrating non-linear systems in each stage. In Bouhamidi and Jbilou,² the effort was made to reduce the computation cost by transforming the linear system into Stein matrix equations, and this linear system arises in Newton’s method. In Cavaglieri and Bewley,³ it was pointed out that implicit/explicit Runge-Kutta schemes were effective for time marching ordinary differential equations having non-stiff and stiff parts. These problems were discretized using implicit Runge-Kutta for the stiff part, which is often linear, and the use of the explicit RK method for the non-stiff part is often non-linear.

In Cavaglieri and Bewley,³ eight implicit/explicit low storage Runge-Kutta schemes were developed and characterized with better stability than the existing only low-storage implicit/explicit Runge-Kutta scheme. Based on the work of Kennedy and Carpenter,¹ several general-purpose, optimal diagonally implicit Runge-Kutta methods have been presented. All methods were L-stable and stiffly accurate. If there could be the possibility, then many were internally L-stable on stages. An embedded method was also included for the facilitation of step size control. In Kennedy and Carpenter,⁴ five new single diagonally implicit RK methods were presented in an explicit first stage. It was claimed that the presented sixth-order stage order two L stable 6(5) pair explicit first stage, single diagonally implicit RK method was the first of its kind. An approach was given in Kutniv et al.⁵ to construct high-order step numerical methods for solving initial value problems on a finite domain. A transformation was given that reduced the problem on semi-infinite interval and then used the finite irregular grid on a semi-infinite interval. Runge-Kutta methods and Taylor series were developed. Special two-derivative Runge-Kutta type methods have been introduced in Lee et al.⁶ for finding the solution of third-order ordinary differential equations, which involved the fourth derivative of the solution. B series theory and rooted tree theory were proposed for deriving order conditions for the special two derivative RK methods. For the effectiveness and accuracy of proposed methods, various test problems were given, and also comparison was made to existing methods. In Rainwater and Tokman,⁷ exponential propagation iterative Runge-Kutta methods were extended to construct split exponential propagation RK integrators for the system of ODEs. It was demonstrated that the derivation of order conditions for the proposed integrators was followed from the bicolor trees-based B series. It was shown in Rainwater and Tokman⁷ that the specific scheme could be custom-built to improve the computational efficiency of the problem. For more numerical methods, readers can see Chen et al.⁸, Figueroa et al.,⁹ Murua,¹⁰ Qin and Zhu,¹¹ Simos,¹² You and Chen,¹³ and references therein.

On the other hand, the literature has extensively discussed flows over a stretching sheet with different kinds of fluid. Not only Newtonian, some non-Newtonian, laminar incompressible fluid flows with heat transfer effects have also been explored in the existing literature. In most cases, the governing equations of the heat and mass transfer in fluid flow have been expressed mathematically in partial differential equations and then further reduced into ordinary differential equations. Some shooting methods, homotopy analysis methods, and Matlab solvers have been utilized to solve the boundary value problems. The effect of magnetic field on forced convective power-law fluid has been investigated in Ahmed and Iqbal.¹⁴ This flow was studied in the annular sector duct. The control volume approach solved the reduced equations, and successive over-relaxation or strongly implicit procedures were adopted to solve resulting algebraic equations. The choice of the method was based on the flow behavior index. Laminar, electrically conducting, incompressible non-Newtonian power-law fluid has been studied in Shamshuddin et al.¹⁵ over an exponentially stretching sheet, and power-law slip velocity was also considered. The spectral quasi-linearization method was adopted to solve the obtained set of non-linear ordinary differential equations. It was found from the plot of obtained results that temperature de-escalated by enhancing Hall current. In Prasad et al.,¹⁶ viscous, incompressible steady state, electrically conducting MHD power-law fluid flow passing through vertically stretching sheet was considered. The transformed set of boundary value problems were solved numerically using the Keller-Box method. The velocity and temperature profile behavior was discussed for parameters involved in the transformed set of non-linear ordinary differential equations. The MHD electrically conducting mixed convective fluid flow past a stretching surface has been studied¹⁷ with heat absorption/generation effects. The set of ordinary differential equations has been tackled by the implicit finite-difference method for coupled non-similar flow. The solution depended on the magnetic parameter, power-law fluid index, modified Richardson number, the generalized Prandtl number, and heat generation and radiation parameters. The analysis of power-law fluid flow over radially stretching sheets has been analyzed¹⁸ with heat transfer effects. The governing unsteady boundary layer equations in partial differential equations have been converted into non-linear ordinary differential equations using new similarity transformations. The system of equations was solved by the homotopy analysis method and shooting method. It was found that the thickness of the boundary layer for velocity and temperature was a decaying function of the unsteadiness parameter. More about power-law fluid, readers can see Abd El-Aziz and Saleem¹⁹, Prasad et al.²⁰, Srinivasacharya and Swamy Reddy²¹, and references therein.

This work aims to present a two-stage numerical scheme for solving ordinary differential equations. The first stage of the scheme is explicit. In contrast, the second stage is implicit, and the scheme’s accuracy is achieved by comparing coefficients in Taylor series expansions. First, the scheme’s stability region is found when applied to the linear differential equation. Later on, the scheme is applied to non-linear boundary value problems obtained using transformations on the governing equations of heat and mass transfer of power-law fluid flow over the sheet. Finally, since equations are non-linear and are linearized by finding the Jacobian of the system and then using the Gauss-Seidel iterative procedure, the whole set of equations are solved and compared with the non-linear approach and existing standard finite difference methods.

For the construction of the numerical scheme, consider the following equation having the form

y' = f (y)

(1)

The first stage of the scheme is expressed as

{\bar{y}}_{i + 1} = y_{i} + {hy}_{i}'

(2)

Where $h$ is the step size. This stage (2) is explicit, and the implicit second stage can be expressed as

y_{i + 1} = y_{i} + h ({ay}_{i + 1}' + b {\bar{y}}_{i + 1}' + {cy}_{i}')

(3)

Stage (3) contains three unknown parameters $a, b,$ and $c$ and values of these three unknown parameters are found by applying the Taylor series expansions. Therefore expanding $y_{i + 1}$ and $y_{i + 1}'$ as

y_{i + 1} = y_{i} + {hy}_{i}' + \frac{h^{2}}{2} y_{i}^{″} + \frac{h^{3}}{6} y_{i}^{″'} + O (h^{4})

(4)

y_{i + 1}' = y_{i}' + {hy}_{i}^{″} + \frac{h^{2}}{2} y_{i}^{″'} + O (h^{3})

(5)

Substituting equation (2) into equation (3) results

y_{i + 1} = y_{i} + h [{ay}_{i + 1}' + b {y_{i}' + {hy}_{i}^{″}} + {cy}_{i}']

(6)

Substituting equations (4) and (5) into (6) gives

\begin{matrix} y_{i} + {hy}_{i}' + \frac{h^{2}}{2} y_{i}^{″} + \frac{h^{3}}{6} y_{i}^{″'} + \dots = y_{i} \\ + h [a {y_{i}' + {hy}_{i}^{″} + \frac{h^{2}}{2} y_{i}^{″'} + \dots} + b {y_{i}' + {hy}_{i}^{″}} + {cy}_{i}'] \end{matrix}

(7)

Comparing coefficients of ${hy}_{i}', h^{2} y_{i}^{″}$ and $h^{3} y_{i}^{″'}$ on both side of equation (7) gives

1 = a + b + c,

(8)

\frac{1}{2} = a + b,

(9)

\frac{1}{6} = \frac{a}{2}

(10)

Equations (8)–(10) are linear, containing three parameters $a, b,$ and $c$ . Solving these equations (8)–(10) gives the values of unknown parameters as

a = \frac{1}{3}, b = \frac{1}{6}, c = \frac{1}{2} .

(11)

Thus, the second stage of the proposed scheme can be expressed as

y_{i + 1} = y_{i} + h (\frac{1}{3} y_{i + 1}' + \frac{1}{6} {\bar{y}}_{i + 1}' + \frac{1}{2} y_{i}')

(12)

This stage (12) of the proposed scheme is implicit, and since the terms of the Taylor series up to third order are compared, the resulting scheme is third-order accurate. Finally, the application of the scheme is given for two examples.

Example 1: Heat and mass transfer of power-law fluid flow

Consider a steady two-dimensional incompressible non-Newtonian power-law laminar fluid flow over the plate. Let the plate be placed in $xz$ -plane, where $x$ - axis is considered to be horizontal, and $y$ -axis is perpendicular to the $x$ -axis. The fluid is placed in space $y > 0$ . The flow moves toward the positive $x$ -axis when the plate is stretched. Let the fluid be electrically conductive and a magnetic field having the strength $B_{0}$ is applied perpendicular to the flow. The stretching velocity of the plate is $U_{w} = ax$ . The geometry of the problem is shown in Figure 1. By considering the effects of thermal radiations and chemical reaction, the governing equations of this phenomenon can be expressed as

\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0

(13)

u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = - ν \frac{\partial}{\partial y} {(\frac{\partial u}{\partial y})}^{n} - \frac{σ B_{0}^{2}}{ρ} u

(14)

u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} = \frac{k}{ρ c_{p}} \frac{\partial^{2} T}{\partial y^{2}} - \frac{1}{ρ c_{p}} \frac{\partial q_{r}}{\partial y}

(15)

u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = D \frac{\partial^{2} C}{\partial y^{2}} - k_{1} (C - C_{\infty})

(16)

Figure 1.

Geometry of the problem.

Subject to the boundary conditions

\begin{matrix} u - U_{w}, v = 0, T = T_{w}, C = C_{w} when y = 0 \\ y = 0, u \to 0, T \to T_{\infty}, C \to C_{\infty} when y \to \infty \end{matrix}}

(17)

In literature, the linearized Roseland flux can be expressed as

q_{r} = - \frac{16 σ^{*} T_{\infty}^{3}}{3 k^{*}} \frac{\partial T}{\partial y}

(18)

Following transformations are considered that can be used to reduce the set of partial differential (13)–(17) into a set of ordinary differential equations

\begin{matrix} u = a x f^{'}, v = - \frac{a}{x (n + 1)} {(\frac{a^{2 - n}}{ν x^{2 n}})}^{- 1 / n + 1} [2 n f - (n - 1) η f^{'}] \\ η = \frac{y}{x} {(\frac{a^{2 - n} x^{2}}{ν})}^{\frac{1}{n + 1}}, θ = \frac{T - T_{\infty}}{T_{w} - T_{\infty}}, ϕ = \frac{C - C_{\infty}}{C_{w} - C_{\infty}} \end{matrix}}

(19)

Substituting transformations (19) into equations (13)–(17) gives the set of ordinary differential equations given as

f^{' 2} - \frac{2 n}{n + 1} f f^{″} = n {(- f^{″})}^{n - 1} f^{″'} - Mf'

(20)

- \frac{2 n}{n + 1} f θ' = \frac{1}{N P_{r}} (1 + \frac{4}{3} R_{d}) θ ″

(21)

- \frac{2 n}{n + 1} f ϕ' = \frac{1}{N S_{c}} ϕ^{″} - γ ϕ

(22)

Subject to the boundary conditions

\begin{matrix} f = 0, f' = 1, θ = 1, ϕ = 1 when η = 0 \\ f' \to 0, θ \to 0, ϕ \to 0 when η \to \infty \end{matrix}}

(23)

Where $M$ is the magnetic parameter, $N P_{r}$ is the modified Prandtl number, $N S_{c}$ is the modified Schmidt number, and $γ$ be dimensionless reaction rate parameter and these are defined as

\begin{matrix} M = \frac{σ B_{°}^{2}}{ρ a}, N P_{r} = \frac{a x^{2}}{α} {(R e_{x})}^{\frac{2}{n + 1}}, N S_{c} = \frac{a x^{2}}{D} {(R e_{x})}^{\frac{2}{n + 1}}, \\ γ = \frac{k}{a}, R_{d} = \frac{4 σ^{*} T_{\infty}^{3}}{k k^{*}} \end{matrix}

where $R e_{x} = \frac{a^{2 - n} x^{n}}{ν}$ .

The quantities of physical interest local Nusselt and Sherwood numbers are defined as

\begin{matrix} N u_{x} = \frac{x q_{w}}{k (T_{w} - T_{\infty})} \\ S h_{x} = \frac{x q_{m}}{D (C_{w} - C_{\infty})} \end{matrix}

where $q_{w}$ and $q_{m}$ are called surface heat and mass fluxes, respectively, and for the considered problem, these are defined as

\begin{matrix} q_{w} = - k {\frac{\partial T}{\partial y} |}_{y = 0} - \frac{16 σ^{*} T_{\infty}^{3}}{3 k^{*}} {\frac{\partial T}{\partial y} |}_{y = 0} \\ q_{m} = - D {\frac{\partial C}{\partial y} |}_{y = 0} \end{matrix}

The non-dimensional forms of local Nusselt and Sherwood numbers are given as

{Re}_{x}^{- \frac{1}{n + 1}} N u_{x} = - (1 + \frac{4}{3} R_{d}) θ' (0)

(24)

{Re}_{x}^{- \frac{1}{n + 1}} S h_{x} = - \emptyset' (0)

(25)

Numerical solution

To solve equations (20)–(23), a proposed scheme of order three is employed. For this reason, the system of second and third-order ODEs are reduced into a system of first-order ODEs in the following manner.

Let

f' = f_{1}

(26)

and

f_{1}' = f_{2}

(27)

It implies

f_{2}' = \frac{1}{n {(- f_{2})}^{n - 1}} {f_{1}^{2} - \frac{2 n}{n + 1} f f_{2} + M f_{1}}

(28)

Also let

θ' = θ_{1}

(29)

θ_{1}' = N P_{r} (\frac{3}{3 + 4 R_{d}}) (- \frac{2 n}{n + 1} f θ_{1})

(30)

Let

ϕ' = ϕ_{1}

(31)

ϕ_{1}' = N S_{c} (- \frac{2 n}{n + 1} f ϕ_{1} + γ ϕ)

(32)

Subject to the given and assumed initial conditions

\begin{matrix} f (0) = 0, f_{1} (0) = 1, f_{2} (0) = x_{1}, \\ θ (0) = 1, θ_{1} (0) = x_{2}, ϕ (0) = 1, ϕ_{1} (0) = x_{3} \end{matrix}

(33)

Where $x_{1}, x_{2}$ , and $x_{3}$ are unknowns, and later on, the values of these unknowns will be found by shooting approach using known boundary conditions.

Since equations (28), (30), and (32) are non-linear ODEs and these are linearized by finding the Jacobian evaluated at the previous grid point. The system of equations (26)–(32) can be expressed as

Y' = F + J Y

(34)

$Y = {[f, f_{1}, f_{2}, θ, θ_{1}, ϕ, ϕ_{1}]}^{t}$ and

F = [\begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} f_{1} \\ f_{2} \\ \frac{1}{n {(- f_{2})}^{n - 1}} {f_{1}^{2} - \frac{2 n}{n + 1} f f_{2} + M f_{1}} \end{matrix} \\ θ_{1} \end{matrix} \\ N P_{r} (\frac{3}{3 + 4 R_{d}}) (- \frac{2 n}{n + 1} f θ_{1}) \end{matrix} \\ ϕ_{1} \end{matrix} \\ N S_{c} (- \frac{2 n}{n + 1} f ϕ_{1} - γ ϕ) \end{matrix}]

and the Jacobian is expressed as

J = [\begin{matrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ - \frac{2 {(- f_{2})}^{1 - n} f_{2}}{n + 1} & \frac{{(- f_{2})}^{1 - n} (2 f_{1} + M)}{n} & J_{3, 3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ \frac{- 2 n θ_{1}}{(n + 1) δ} & 0 & 0 & 0 & \frac{- 2 fn}{(n + 1) δ} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \frac{- 2 nN S_{c} ϕ_{1}}{(n + 1)} & 0 & 0 & 0 & 0 & - γ N S_{c} & \frac{- 2 nfN S_{C}}{(n + 1)} \end{matrix}]

Where $J_{3, 3} = - \frac{2 {(- f_{2})}^{1 - n} f}{n + 1} - \frac{{(- f_{2})}^{- n} (1 - n) (f_{1}^{2} + f_{1} M - \frac{2 f f_{2} n}{1 + n})}{n}$ , and $δ = \frac{1}{N P_{r}} (1 + \frac{4}{3} R_{d})$ .

Applying the first stage of the proposed scheme on (34) yields

{\bar{Y}}_{i + 1} = Y_{i} + h (F_{i} + J_{i} Y_{i})

(35)

and the second stage of the scheme is expressed as

\begin{matrix} Y_{i + 1} = Y_{i} + h \\ [\frac{1}{3} (F_{i} + J_{i} Y_{i + 1}) + \frac{1}{6} (F_{i} + J_{i} {\bar{Y}}_{i + 1}) + \frac{1}{2} (F_{i} + J_{i} Y_{i})] \end{matrix}

(36)

For $j = 1, 2, \dots, 7$ and $1 \leq M \leq 7$ , the $m$ th component in equation (36) can be expressed as

\begin{matrix} Y_{m, i + 1} = Y_{m, i} + h F_{m, i} + h \\ [\frac{1}{3} \sum_{j = 1}^{m} J_{(m, j), i} Y_{j, i} + \frac{1}{3} \sum_{j = m + 1}^{7} J_{(m, j), i} Y_{j, i} + \frac{1}{6} \sum_{j = 1}^{7} J_{(m, j), i} {\bar{Y}}_{j, i + 1} + \frac{1}{2} \sum_{j = 1}^{7} J_{(m, j), i} Y_{j, i}] \end{matrix}

(37)

where $J_{(m, j), i}$ represents $j$ th entry of $m$ th line of Jacobian evaluated at grid point $i$ . Employing Gauss-Seidel iterative procedure on equation (37) results,

\begin{matrix} Y_{m, i + 1}^{k + 1} = Y_{m, i}^{k + 1} + {hF}_{m, i}^{k + 1} + h \\ [\frac{1}{3} \sum_{j = 1}^{m} J_{(m, j), i}^{k + 1} Y_{j, i + 1}^{k + 1} + \frac{1}{3} \sum_{j = m + 1}^{7} J_{(m, j), i}^{k + 1} Y_{j, i + 1}^{k} + \frac{1}{6} \sum_{j = 1}^{7} J_{(m, j), i}^{k + 1} {\bar{Y}}_{j, i + 1}^{k + 1} + \frac{1}{2} \sum_{j = 1}^{7} J_{(m, j), i}^{k + 1} Y_{j, i}^{k + 1}] \end{matrix}

(38)

and equation (38) can be expressed as

{\bar{Y}}_{m, i + 1} = Y_{m, i}^{k + 1} + h [F_{m, i}^{k + 1} + \sum_{j = 1}^{7} J_{(m, j), i}^{k + 1} Y_{j, i}^{k + 1}]

(39)

Stability analysis

The stability analysis of linear differential equations for the scalar case can be found in Appendix. In this section, the stability of the homogeneous equation (34) is given. The homogeneous equation is expressed as

Y' = J Y

(40)

Applying the first stage of the proposed scheme using the Gauss-Seidel iterative method yields equation (39) with $F_{m, i}^{k + 1} = 0, \forall m$ . Also, the second stage of the proposed scheme for equation (40) results from equation (38) with $F_{m, i}^{k + 1} = 0, \forall m$ . Here, it is assumed that all components of the system behave similarly and so,

\begin{matrix} Y_{m, i + 1}^{k + 1} = E^{k + 1} e^{(i + 1) I ψ}, Y_{m, i}^{k + 1} = E^{k + 1} e^{iI ψ} \\ Y_{m, i + 1}^{k} = E^{k} e^{(i + 1) I ψ} \\ {\bar{Y}}_{m, i + 1} = {\bar{E}}^{k + 1} e^{(i + 1) I ψ} for 1 \leq m \leq 7 \end{matrix}}

(41)

Where $I = \sqrt{- 1}$ . Substituting transformation (41) into equation (39) with $F_{m, i}^{k + 1} = 0$ yields

\begin{matrix} {\bar{E}}^{k + 1} e^{(i + 1) I ψ} = E^{k} e^{(i + 1) I ψ} + h [\sum_{j = 1}^{7} J_{(m, j), i}^{k + 1} E^{k + 1} e^{iI ψ}] \\ = E^{k + 1} e^{iI ψ} + h E^{k + 1} e^{iI ϕ} \sum_{j = 1}^{7} J_{(m, j), i}^{k + 1} \end{matrix}

(42)

Substituting transformations (41) into equation (38) with $F_{m, i}^{k + 1} = 0 \forall m$ results

\begin{matrix} (1 - \bar{a}) E^{k + 1} e^{(i + 1) I ψ} = E^{k + 1} e^{iI ψ} + h \\ [\frac{1}{3} E^{k + 1} e^{(i + 1) I ψ} \bar{b} + \frac{1}{3} E^{k} e^{(i + 1) I ψ} \bar{c} + \frac{1}{6} \bar{E} e^{(i + 1) I ψ} \bar{d} + \frac{1}{2} E^{k + 1} e^{iI ψ} \bar{d}] \end{matrix}

(43)

Where $\bar{a} = \frac{h}{3} J_{(m, m), i}^{k + 1}, \bar{b} = \sum_{j = 1}^{m - 1} J_{(m, j), i}^{k + 1}, \bar{c} = \sum_{j = m + 1}^{7} J_{(m, j), i}^{k + 1}, \bar{d} = \sum_{j = 1}^{7} J_{(m, j), i}^{k + 1}$

Substituting equation (42) into equation (43) and dividing both sides of the resulting equation by $e^{iI ψ}$ gives

\begin{matrix} (1 - \bar{a}) E^{k + 1} e^{I ψ} = E^{k + 1} + h \\ [\frac{1}{3} E^{k + 1} e^{I ψ} \bar{b} + \frac{1}{3} E^{k} e^{I ψ} \bar{c} + \frac{1}{6} (E^{k + 1} + h E^{k + 1} \bar{d}) + \frac{1}{2} E^{k + 1} \bar{d}] \\ [(1 - \bar{a}) e^{I ψ} - 1 - \frac{h \bar{b}}{3} e^{I ψ} - \frac{h}{6} - \frac{h^{2}}{6} \bar{d} - \frac{h}{2} \bar{d}] E^{k + 1} = \frac{\bar{c}}{3} e^{I ψ} E^{k} \end{matrix}

And the amplification factor is expressed as

\frac{E^{k + 1}}{E^{k}} = \frac{\bar{c} / 3 e^{I ψ}}{(1 - \bar{a}) e^{I ψ} - 1 - \frac{h \bar{b}}{3 e^{I ψ}} - \frac{h}{6} - \frac{h^{2}}{6} \bar{d} - \frac{h}{2} \bar{d}}

The stability condition can be expressed as

| \frac{\bar{c} / 3 e^{I ψ}}{(1 - \bar{a}) e^{I ψ} - 1 - \frac{h \bar{b}}{3} e^{I ψ} - \frac{h}{6} - \frac{h^{2}}{6} \bar{d} - \frac{h}{2} \bar{d}} | < 1

This implies

| \bar{c} / 3 e^{I ψ} | < | (1 - \bar{a}) e^{I ψ} - 1 - \frac{h \bar{b}}{3} e^{I ψ} - \frac{h}{6} - \frac{h^{2}}{6} \bar{d} - \frac{h}{2} \bar{d} |

Using different step sizes by fixing $\bar{a}, \bar{b}, \bar{c}$ , and $\bar{d}$ with $ψ = ψ_{s} = \frac{s π}{M - 1} \forall s = 0, 1, 2, \dots M - 1$ where $M$ represents the number of grid points. The stability conditions can be determined for particular cases. For detail of this type of stability, readers can see Mazumder.²²

Results and discussions

The proposed numerical scheme has been applied to the linear ordinary differential equation given in the first example. Also, it is applied for solving the problem of heat and mass transfer of boundary layer flow. The proposed numerical scheme for the non-linear problem has been applied in two different directions. In one direction, the system of non-linear differential equations is linearized by finding the Jacobian of the whole system of first-order differential equations. The Jacobian is evaluated at the previous grid point. Then, the Gauss-Seidel iterative method is employed to solve the difference equation obtained by applying the proposed scheme. Because the proposed strategy applies to first-order differential equations, the shooting method is considered. Therefore, the scheme is based on three different approaches. The proposed scheme performs the discretization, then the Gauss-Seidel iterative method is considered, and for finding missing initial conditions, a shooting approach is considered. This whole approach requires two sets of initial guesses to get the solution of the required differential equations.

The second approach is non-linearized instead of linearized. In this non-linearized approach, required differential equations are discretized using the proposed scheme, and a modified iterative method is employed to solve the resulting non-linear difference equations. This approach is more accurate than the linearized one, which will be demonstrated in the residuals plots. Both discretizations and iterative methods converged to the accurate solution for the case of Newtonian fluid with $n = 1$ . Figure 2 show the comparison of linearized (Figure 3, left) and non-linearized (Figure 3, right) approaches by finding residuals over the iterations. This Figure 2 shows the effectiveness of the non-linear approach. All three residuals for three non-linear differential equations are approaching zero approximately when iterations are increased. But all the residuals in the linearized approach do not tend to zero by increasing the iterations. The residuals are found using the procedure given in Nawaz and Shoaib.²³

Figure 2.

Comparison of proposed scheme using iterative procedures for the linearized and non-linearized system using $n = 1, M = 0.7, R_{d} = 0.1, N P_{r} = 2.7, N S_{c} = 1.7, γ = 0.1, N = 90, L = 4$ .

Figure 3.

Comparison of proposed scheme using iterative procedures for the linearized and non-linearized system over independent variable $η$ using $n = 1, M = 0.7, R_{d} = 0.1, N P_{r} = 2.7, N S_{c} = 1.7, γ = 0.1, N = 90, L = 4$ .

In this procedure, one extra numerical derivative is found, and the final residual is computed. The maximum number of iterations consumed by linearized and non-linearized approaches is the same. Figure 3 shows the residuals plots over the independent variable $η$ . This Figure 3 elucidates the comparison of linearized and non-linearized approaches. By looking at this Figure 3, it can be concluded that the non-linearized approach did better than the linearized one. This is the main advantage of the non-linear approach. But the linearized approach provides the facility to perform some stability analysis. Figures 4 and 5 show the comparisons of the proposed scheme with the classical/standard second-order finite difference methods. The classical standard finite difference method with the trapezoidal//Keller-Box method is employed to discretize the non-linear ordinary differential equations. Afterward, an iterative method is employed to solve the resulting difference equations. The third-order non-linear ordinary differential equation is reduced to a system of first and second-order ordinary differential equations. The trapezoidal/Keller-Box method is employed on first-order differential equations, and standard/classical second-order central difference formulas are employed to discretize second-order differential equations. From this comparison given in Figures 4 and 5, it can be observed that the proposed scheme performs better than the classical finite difference method. The observation can be made from both types of comparisons, the one made on the consumption of the maximum number of iterations, and the other is made on the smaller maximum residual. The huge difference between the consumption of iterations from existing and proposed schemes can be seen. Figures 6 to 8 show respectively velocity, temperature, and concentration profiles with the variation of parameter $n$ . The behavior of the velocity profile is dual, whereas the temperature and concentration profiles de-escalate by enhancing parameter $n$ . Figure 9 elucidates the impact of magnetic parameters on the velocity profile. Velocity de-escalates with rising values of magnetic parameter $M$ . Since applying magnetic force in the electrically conducting fluid generates Lorentz force, this force resists the fluid’s velocity, which results in a de-escalation in the velocity profile. Figure 10 deliberates the impact of radiation parameter $R_{d}$ on the temperature profile. Temperature escalates by enhancing the values of radiation parameters because radiations in the fluid enhance the amount of energy, which is responsible for increasing the fluid’s temperature. Figure 11 elucidates the impact of modified Prandtl number $N P_{r}$ on temperature profile. Temperature decays with rising values of modified Prandtl number $N P_{r}$ because growth in Prandtl number produces decay in thermal conductivity, which is responsible for de-escalating the thermal conductivity, so temperature decreases, the impact of modified Schmidt number on concentration profile is depicted in Figure 12. The concentration profile de-escalates by enhancing the modified Schmidt number. It happens because an increase in modified Schmidt number produces decay in mass diffusivity. This decay means less deriving force for diffusion applied on concentration by considering Fick’s law. The impact of reaction rate parameter $γ$ on concentration profile can be seen in Figure 13. From this Figure 13, it can be observed that the concentration profile de-escalates by increasing the reaction rate parameter. Since the involvement of chemical reaction results from the reaction in concentration profile that leads to breaking of chemical bonds in atoms which slows down the concentration profile. In most of the Figures 2 to 13, $N$ denotes the maximum number of grid points, and $L$ denotes the length of the domain. The ranges of considered parameters depend on the physical situations but numerically, the ranges of parameters depend on the convergence of the constructed Matlab code. Numerically any value of the parameter can be considered on which the code satisfies the given criteria of convergence. So, the code will be stopped if norms of differences for solutions of each equation computed on two consecutive iterations are less than the value near zero. After convergence of the code, plots have been drawn. The stopping criteria for each equation can be expressed as

‖ g^{k + 1} - g^{k} ‖_{2} < ϵ

Figure 4.

Comparison of proposed and existing numerical methods over iterations using $n = 1, M = 0.7, R_{d} = 0.1, N P_{r} = 2.7, N S_{c} = 1.7, γ = 0, N = 90, L = 4$ .

Figure 5.

Comparison of proposed and existing numerical methods over independent variable $η$ using $n = 1, M = 0.7, R_{d} = 0.1, N P_{r} = 2.7, N S_{c} = 1.7, γ = 0, N = 90, L = 4$ .

Figure 6.

Velocity profile with the variation of $n$ using $M = 0, N = 60$ .

Figure 7.

Temperature profile with the variation of $n$ using $N = 60, M = 0.1, R_{d} = 0, N P_{r} = 1$ .

Figure 8.

Concentration profile with the variation of $n$ using $N = 60, L = 7, M = 0.4, N S_{c} = 5, γ = 0.7$ .

Figure 9.

Velocity profile with the variation of magnetic parameter $M$ using $n = 1, N = 90$ .

Figure 10.

Temperature profile with the variation of radiation parameter $R_{d}$ using $n = 1, N = 90, M = 0.5, N P_{r} = 1$ .

Figure 11.

Temperature profile with the variation of modified Prandtl number $N P_{r}$ using $n = 1, N = 90, M = 0.5, R_{d} = 0.5$ .

Figure 12.

Concentration profile with the variation of modified Schmidt number $N S_{c}$ using $n = 1, N = 90, M = 0.5, γ = 0.9$ .

Figure 13.

Concentration profile with the variation of reaction rate parameter $γ$ using $n = 1, N = 90, M = 0.5, N S_{c} = 1.7$ .

Where $g$ is a vector of solution for any equation with dependent variables $f$ , $θ$ , and $\emptyset$ , $k$ represents the iteration number, and $ϵ$ denotes the positive value near zero.

Table 1 shows numerical values for local Nusselt and local Sherwood numbers with different contained parameter values variations. For solving the considered model of non-Newtonian power-law fluid, the Matlab code found complex values with small imaginary parts. But code uses real parts of complex entries and ignored the imaginary parts. These complex values are found for the non-Newtonian case when $n \neq 1,$ but for Newtonian fluid, real values are found, but due to the use of shooting technique, an updated guess is provided to Matlab built-in solver fsolve.

Table 1.

List of numerical values for local Nusselt and local Sherwood numbers.

$n$	$M$	$NPr$	$Rd$	$NSc$	$γ$	${Re}_{x}^{- \frac{1}{n + 1}} N u_{x}$	${Re}_{x}^{- \frac{1}{n + 1}} S h_{x}$
$0.3$	$0.1$	$1$	$0.1$	$3$	$0.1$	$0.3648$ ^*	$0.8934$ ^*
$0.7$						$0.5395$ ^*	$1.1652$ ^*
$1$						$0.5962$	$1.2739$
	0.4					$0.5679$	$1.2498$
		$4$				$1.4046$	$1.2498$
			$0.4$			$1.5797$	$1.2498$
				$5$		$1.5797$	$1.6733$
					$0.7$	$1.5797$	$2.3962$

The Matlab code ignored imaginary values.

Figures 14 to 18 are drawn using software used for solving partial differential equations, and it uses the finite element method to solve the required equations. This software is used for solving a model of heat and mass transfer of power-law fluid with effects of reaction in concentration. The variable viscosity for power-law fluid is used as

\begin{matrix} μ = m {({\overset{\cdot}{γ}}_{1})}^{n - 1} \\ \dot{{\overset{\cdot}{γ}}_{1} = \max (\sqrt{D : D}, {\overset{\cdot}{γ}}_{1, \min})}, D = \frac{1}{2} [\nabla u + {(u)}^{T}] \end{matrix}

Figure 14.

(a) Velocity surface, $(b)$ streamlines, $(c)$ temperature surface, $(d)$ isothermal contours using $u_{w} = 0.01, m = 0.01, n = 1.5$ .

Figure 15.

(a) Velocity surface, $(b)$ streamlines, $(c)$ temperature surface, $(d)$ isothermal contours using $u_{w} = 0.01, m = 0.01, n = 0.5$ .

Figure 16.

(a) Velocity surface, $(b)$ streamlines, $(c)$ temperature surface, $(d)$ isothermal contours using $u_{w} = 0.1, m = 0.01, n = 0.5$ .

Figure 17.

(a) Velocity surface, $(b)$ streamlines, $(c)$ temperature surface, $(d)$ isothermal contours using $u_{w} = 0.4, m = 0.01, n = 0.5$ .

Figure 18.

Concentration surfaces using $u_{w} = 0.1, m = 0.01, n = 2 (a) d_{c} = 1 (b) d_{c} = 0.01 (c) d_{c} = 0.0001 (d) d_{c} = 0.0000001$ .

The impact of flow behavior index $n$ on velocity and temperature profiles can be seen in Figures 14 to 15. Velocity and temperature escalate by enhancing the value of the parameter $n$ . The impact of wall velocity $u_{w}$ on velocity and temperature profiles can be seen in Figures 16 to 17. From these Figures 16 and 17, it can be observed that the velocity initially increases by enhancing the wall velocity and leads to growth in the temperature. The impact of mass diffusivity on concentration profile can be visualized in Figure 18. The thickness of the boundary layer grows by the rise of mass diffusivity. The reaction coefficient used in Figure 18 is $R_{c} = - 2$ .

Conclusion

A numerical scheme has been given that can be used to solve first-order linear and non-linear ordinary differential equations. The scheme is third-order accurate, and by finding the stability function, a stability region is plotted using software Mathematica. The mathematical model in the form of partial differential equations for heat and mass transfer of chemically reactive MHD non-Newtonian power-law fluid flow over a stretching sheet with effects of radiations has been presented. Later on, these partial differential equations were converted into a set of non-linear ordinary differential equations using transformations. The non-linear system of equations was linearized by determining the Jacobian. The numerical approach was used to solve both linearized and non-linearized equations using a combination of iterative methods. The comparison showed that the presented numerical scheme was more effective and computationally less expensive than the standard/classical second-order central difference method. The presented numerical scheme can be applied further to solve different ordinary differential equations in science and engineering. Following the completion of this work, it is possible to propose additional applications for the current methodology.^24–28 Additionally, the developed method is simple to implement and may be utilized to solve a broader class of differential equations encountered in practice and theory.

Footnotes

Appendix

Example 2: Consider the linear differential equation of the form

(A1)

y' = ky

Subject to the initial condition

(A2)

y (0) = 1

The exact solution of equation (A1) is $y = e^{kt}$ . Applying the first stage of the scheme to equation (A1) yields

(A3)

{\bar{y}}_{i + 1} = y_{i} + h (k y_{i})

The second stage of the proposed scheme is expressed as

(A4)

y_{i + 1} = y_{i} + h [\frac{1}{3} y_{i + 1}' + \frac{1}{6} {\bar{y}}_{i + 1}' + \frac{1}{2} y_{i}']

Acknowledgements

The authors are thankful to the anonymous referee for their suggestions and helpful comments that improved this article. We are also grateful to Vice-Chancellor, Air University, Islamabad, for providing an excellent research environment and facilities. The second author also thanks Prince Sultan University for funding this research work

Handling Editor: Chenhui Liang

Author contributions

All the authors contributed equally to this research work.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was funded by Prince Sultan University, grant number RG-DES-2017-01-17.

ORCID iD

Muhamamd Shoaib Arif

Code availability

A computational code for the proposed discrete model scheme may be made available to readers upon request for convenience.

Data availability statement

The manuscript included all required data and information for its implementation.

References

Kennedy

Carpenter

MH.

Higher-order additive Runge–Kutta schemes for ordinary differential equations. Appl Numer Math 2019; 136: 183–205.

Bouhamidi

Jbilou

Stein implicit Runge–Kutta methods with high stage order for large-scale ordinary differential equations. Appl Numer Math 2011; 61: 149–159.

Cavaglieri

Bewley

Low-storage implicit/explicit Runge–Kutta schemes for the simulation of stiff high-dimensional ODE systems. J Comput Phys 2015; 286: 172–193.

Kennedy

Carpenter

MH.

Diagonally implicit Runge–Kutta methods for stiff ODEs. Appl Numer Math 2019; 146: 221–244.

Kutniv

Datsko

Kunynets

, et al. A new approach to constructing of explicit one-step methods of high order for singular initial value problems for non-linear ordinary differential equations. Appl Numer Math 2020; 148: 140–151.

Lee

Senu

Ahmadian

, et al. Numerical study of third-order ordinary differential equations using a new class of two derivative Runge-Kutta type methods. Alex Eng J 2020; 59: 2449–2467.

Rainwater

Tokman

A new class of split exponential propagation iterative methods of Runge–Kutta type (sEPIRK) for semilinear systems of ODEs. J Comput Phys 2014; 269: 40–60.

Chen

Shi

Liu

, et al. Trigonometrically fitted two-derivative Runge-Kutta-Nyström methods for second-order oscillatory differential equations. Appl Numer Math 2019; 142: 171–189.

Figueroa

Jackiewicz

Löhner

Efficient two-step Runge-Kutta methods for fluid dynamics simulations. Appl Numer Math 2021; 159: 1–20.

10.

Murua

Runge-Kutta-Nyström methods for general second order ODEs with application to multi-body systems. Appl Numer Math 1998; 28: 387–399.

11.

Qin

Zhu

WJ.

A note on stability of a three-stage difference scheme for ordinary differential equations. Comput Math Appl 1993; 25: 35–44.

12.

Simos

TE.

Modified Runge-Kutta methods for the numerical solution of ODEs with oscillating solutions. Appl Math Comput 1997; 84: 131–143.

13.

You

Chen

Direct integrators of Runge–Kutta type for special third-order ordinary differential equations. Appl Numer Math 2013; 74: 128–150.

14.

Ahmed

Iqbal

MHD power law fluid flow and heat transfer analysis through Darcy Brinkman porous media in annular sector. Int J Mech Sci 2017; 130: 508–517.

15.

Shamshuddin

Salawu

Ogunseye

, et al. Dissipative power-law fluid flow using spectral quasi linearization method over an exponentially stretchable surface with Hall current and power-law slip velocity. Int Commun Heat Mass Transf 2020; 119: 104933.

16.

Prasad

Datti

Vajravelu

Hydromagnetic flow and heat transfer of a non-Newtonian power law fluid over a vertical stretching sheet. Int J Heat Mass Transf 2010; 53: 879–888.

17.

Chen

CH.

Magneto-hydrodynamic mixed convection of a power-law fluid past a stretching surface in the presence of thermal radiation and internal heat generation/absorption. Int J Non Linear Mech 2009; 44: 596–603.

18.

Ahmed

Mahmood

Iqbal

, et al. Axisymmetric flow and heat transfer over an unsteady stretching sheet in power law fluid. J Mol Liq 2016; 221: 386–393.

19.

Abd El-Aziz

Saleem

. Numerical simulation of entropy generation for power-law liquid flow over a permeable exponential stretched surface with variable heat source and heat flux. Entropy 2019; 21: 484.

20.

Prasad

Pal

Datti

PS.

MHD power-law fluid flow and heat transfer over a non-isothermal stretching sheet. Commun Nonlinear Sci Numer Simul 2009; 14: 2178–2189.

21.

Srinivasacharya

Swamy Reddy

Chemical reaction and radiation effects on mixed convection heat and mass transfer over a vertical plate in power-law fluid saturated porous medium. J Egypt Math Soc 2016; 24: 108–115.

22.

Mazumder

. Numerical methods for partial differential equations. The Ohio state university, USA: Academic Press, 2015.

23.

Nawaz

Shoaib Arif

Keller-Box shooting method and its application to nanofluid flow over convectively heated sheet with stability and convergence. Numer Heat Transf B Fundam 2019; 76: 152–180.

24.

Nawaz

Arif

MS.

Modified class of explicit and enhanced stability region schemes: application to mixed convection flow in a square cavity with a convective wall. Int J Numer Methods Fluids 2021; 93: 1759–1787.

25.

Nawaz

Arif

A new class of A-stable numerical techniques for odes: application to boundary layer flow. Therm Sci 2020; 25: 97–1675.

26.

Nawaz

Arif

Nazeer

A class of second-order schemes with application to chemically reactive radiative natural convection flow in a rectangular enclosure. Int J Numer Methods Fluids 2021; 93: 3161–3186.

27.

Nawaz

Shoaib Arif

Bibi

, et al. A finite difference method and effective modification of gradient descent optimization algorithm for MHD fluid flow over a linearly stretching surface. Comput Mater Continua 2020; 62: 657–677.

28.

Nawaz

Arif

Shatanawi

, et al. An explicit fourth-order compact numerical scheme for heat transfer of boundary layer flow. Energies 2021; 14: 3396.