As we know, fluid dynamics is central to studying various other disciplines, such as environmental science, engineering, etc. Also, understanding complex real-world problems demands an insight into the specific characteristics of non-Newtonian fluids. They are the fluids that do not follow the linear relation between the shear stress and the velocity gradient. Casson fluids are a family of viscoplastic materials having yield stress, which is used to model the behaviour of many industrial and biological systems.

Several factors of the fluid, such as temperature, mass diffusivity, etc., affect the behaviour of Casson fluids internally and externally. These factors do not allow an easier way of modelling and recreating the movement of Casson fluids. In addition, Casson fluids frequently flow in systems with spatially varying thermal conductivity and mass diffusivity, necessitating the development of computational techniques that can accurately reflect these fluctuations.

Applications in physical chemistry, metrology, biology, oceanography, astrophysics, plasma physics, etc., all highlight the importance of heat transmission. Liquid distillation, heat exchangers, atomic controller refrigeration, and other technological advances rely heavily on heat transmission. In fluid mechanics, researchers have observed that a certain proportion of mechanical energy is converted into thermal or heat energy due to the resistance generated by viscosity between adjacent fluid layers during their motion. We refer to this as a “switch in internal energy.” First, in his essay [1], Brinkman studied the impact of an internal energy change on capillary flow. Using the impacts of internal energy change and heat transport, Jambal et al. [2] established a power law model and estimated the answer using the finite difference approach. The utilization of nanofluids to improve heat transfer has garnered significant interest among academics in recent years due to its extensive applicability in various industries, such as photonics, electronics, energy production, and transportation [3]. In general, metallic fluids tend to exhibit higher thermal conductivity when compared to non-metallic fluids. Hence, it can be observed that the thermal performance of simple fluids is relatively worse when compared to the thermal performance of metallic nano-sized solid particles dispersed in typical fluids. Nanofluids are formed by introducing microstructural particles into ordinary fluids. These particles, typically composed of metals, carbides, carbon nanotubes, or oxides, have dimensions on the nanometer scale [4, 5]. The nanoparticle composition is crucial in hybrid nanofluids, particularly in enhancing distinctive features such as thermal conductivity. Aziz [6] used the shot method to solve the governing equations, demonstrating the impact of viscous dissipation on an energy equation and the effect of altering thickness on momentum equations.

Nanofluids are the subject of many studies because of their superior conduction qualities that can be achieved through various nanofluid compositions [7, 8]. Herein, we list a few studies conducted along these lines. Nasrin and Alim [9] conducted a numerical study of the heat transmission rate for nanofluids containing dual particles.

Furthermore, a method for simulating micro- and nano-scale fluids has been investigated by Nie et al. [10]. In [11], the writers delve into the theoretical framework of hybrid nanofluids’ heat conduction. The effect of hybrid nanofluids on forced convective heat transfer was estimated statistically by Labib et al. [12].

Investigating fluid flow induced by a horizontally translating surface and its impact on thermo-physical characteristics, such as mass diffusivity, thermal conductivity, and viscosity, is a very captivating subject matter for researchers and scholars. Many studies do not account for or assume that a malleable surface’s thermophysical parameters like conductivity, diffusivity, and viscosity are constant. However, the findings of the experiments show that these thermo-physical properties depend on temperature and concentration, especially in the case of a very large temperature differential. As a result, much attention has been focused on how different thermo-physical factors affect surface stretching. The effect of radiation and thermo-physical factors on the flow of a viscous fluid towards a non-uniform permeable medium was explored by Elbarbary et al. [13]. Saleem studied the effects of various fluid properties on viscous fluid flow through a stretchable medium [14].

Hashim et al. [15] proposed the Willaimson fluid model incorporating nanoparticles, where thermophysical parameters were treated as independent variables. Malik et al. [16] investigated how different fluid properties affected the boundary layer flow of a viscous fluid induced by an expanded cylinder. By assuming exponential functions of temperature for both viscosity and thermal conductivity, Mohiuddin et al. [17] can define the behaviour of a viscoelastic fluid. Second-order fluid flow via a mobile medium in the presence of a heat source/sink was studied by Akinbobola et al. [18], who examined the effect of temperature-dependent physical features of the fluid. Muthucumaraswamy [19, 20] solved the constitutive equations for the flow of a viscous fluid across a non-uniform plate using the Laplace method and variable diffusivity. The 1D -diffusion-advection equation was studied by Jia et al. [21] in two different scenarios: (i) when the thermo-physical characteristics are fixed but the flow velocity is not, and (ii) when the flow velocity and the parameters of the fluid are both changeable. The model was solved, and the resulting outcomes for two scenarios were compared.

The thermo-physical parameters that change with temperature and concentration were studied by Li et al. [22] using the finite difference approach to examine their influence on nonlinear transient responses. The effects of temperature and concentration on the transmission of heat and mass in a viscoelastic fluid flow were examined by Qureshi et al. [23]. Researchers in [24] analyzed Maxwell’s fluid flow model for nanoparticles over a heterogeneous medium, considering thermal effects. Near a vertically moving surface, boundary layer flow is due to cooling and heating impulses [25]. The boundary layer flow around an isothermal, free-moving needle was discussed in [26]. [27] examined the heat transfer parameters of forced convection flow over a non-isothermal thin needle. The Boungirono model of nanofluid flow over a rotating needle was analyzed in [28]. Solving the governing equations involved shooting and fourth-order RK methods. The role of heat production and thermal radiation in MHD The effects of an infinite horizontal sheet on the flow of a Casson fluid in two dimensions were studied in [29].

Animasaun [30] examined how vertically uneven surface roughness affected an unstable mixed convection flow. To investigate the flow’s reaction to a chemical reaction and radiation, he applied the shooting method and quadratic interpolation to the model and then solved it. Shah et al. [32] investigated the effect of the Grashof number on the mixed convection flow of different fluids travelling along different surfaces when heat generation was present. In their study, Animasaun et al. [33] looked at how a chemical process, including quartic autocatalysis, might alter the trajectories of various airborne dust particles. Runge-Kutta, shooting, and bvp4c were used to solve the model’s constitutive equations. The influence of nanoparticles’ random mobility in three-dimensional flow was investigated in a recent meta-analysis by Animasaun et al. [34]. They calculated the heat transfer rate due to the Brownian motion of nanoparticles by considering radiation from the surface and local and mass convection. Researchers at [35] examined how alumina nanoparticles behaved in three dimensions when they carried water or were subject to Lorentz force. They looked into how various dimensionless factors affected the velocities involved.

The article [36] delves into heat transfer in Jeffery-Hamel hybrid nanofluid flows involving non-parallel plates. Molybdenum disulfide nanoparticles are suspended in fluids subjected to magnetic fields, heat radiation, and viscous dissipation. The researcher studied the flow of micropolar fluids across a vertical Riga sheet [37]. We look at the nonlinear stretching sheet. A magnetohydrodynamic (MHD) pair stress hybrid nanofluid on a contracting surface is studied in terms of its radiative properties and overall stability [38].

The difficulties in simulating Casson fluids with non-constant thermal and mass diffusivities can be mitigated with the help of stochastic numerical techniques. To capture the inherent stochasticity in real-world systems, these methods add probabilistic features to account for uncertainties and fluctuations in material qualities. Randomness can be due to material contamination, temperature difference or mass concentration change.

In the past, problems with intricate flows were analyzed by finite difference or finite element methods or by the CFD (Computational Fluid Dynamics) simulations to obtain a better viewpoint of fluid dynamics. Although these methods have enhanced our understanding of the subject, they have often failed to reproduce the inherently stochastic behaviours found in numerous real-world systems accurately. The natural uncertainties in several physical processes in fluid flow are due to numerous boundary conditions, material qualities, and environmental effects. If we do not consider these random variables, then the resulting description of the phenomenon may lead to a misleading picture of reality.

In this work, we examine and assess stochastic numerical methodology for modelling of dynamics of the Casson fluid with arbitrary temperature and density gradients for a better view of how uncertainties and fluctuations in material qualities influence the flow of Casson fluid; we will include a stochastic ingredient in the numerical simulations.

Applications in chemical engineering, geophysics, and biomedicine can significantly profit from gaining exact forecasts of the fundamental behaviour of the fluids to optimize the given processes, develop equipment or know the system’s biological behaviour.

We are at the dawn of applying stochastic probability in fluid mechanics; there is a long way to go. The present article goes into this notion. Let’s consider using stochastic predicting in computational fluid mechanics. We will understand mathematical predicting to describe the behaviours of a physical system’s system within which it operates. Computational models need optimization, design, and updating due to external effects like fluctuations in the natural system and internal elements like uncertainty in the model itself.

Numerous scholars are working hard to figure out stochastic partial differential equations and their numerical solutions. Tessitoe [39] made a seminal discovery in this area when he found that linear and infinite-dimensional stochastic differential equations satisfy the same general conditions as the modified solution. The authors of [40] examined the classical form of the stochastic equation under the assumption of homogeneous Dirichlet boundary conditions. The group set out to see if there were any non-trivial positive global solutions and whether or not those solutions were likely to explode in finite time. Researchers in Ref. [41] examined the Holder continuous coefficient obtained with constant coloured noise to study the stochastic partial differential equation (SPDE). Solving a backward double stochastic differential equation (SDE) allows for path-wise uniqueness and deterministic manipulation of the Laplacian. The solution to a system of stochastic differential equations (SDEs) is found by taking weak limits of a sequence of variables. We obtain this sequence by substituting the discrete Laplacian operator for the random variable in the stochastic partial differential equation (SPDE). Altmeyer et al. explained cellular repolarization using a stochastic variant of the Meinhardt equation. The driving noise process has been shown to influence the evolution of solution patterns for stochastic partial differential equations (SPDEs), and such solutions exist [42]. The solution is fully described in the works mentioned above.

Numerical estimation of stochastic partial differential equations (SPDEs) is a formidable challenge. Instead, Gyorgy et al. [43] worked to construct lattice approximations for elliptic stochastic partial differential equations (SPDEs). For white noise on a restricted domain in R d , d = 1 , 2 , 3 , the convergence rate of approximations is calculated. In [44], we look at how to approximate answers to stochastic partial differential equations of the Itô type. The consistency and stability of these approximations with respect to their mean-square error are established by employing explicit and implicit finite difference techniques. The stochastic Fitz-Hugh-Nagumo model was defined, and a numerical solution was given in [45]. This examination shows how well the technique holds up in a Von Neumann environment [46]. investigated the reliability and robustness of the forward Euler method for evaluating stochastic nonlinear advection-diffusion models. In [47], they consider white noise’s spectral power distribution functions and estimate the numerical approximations for the linear, elliptic, and parabolic cases. The approximations of these cases are evaluated using the finite element and difference methods. The relevant literature dealt with the integral approximation techniques, the finite element methods in these contexts, and the weak SPDE formulation.

This research paper proposes a new and novel numerical scheme for solving the problems of unstable non-Newtonian mixed convection flow of fluid with heat and mass transport with the effect of temperature and concentration fluctuations. The proposed methodology combines stochastic methods in a finite difference scheme, which enables the capture of the random behaviour of the fluid flow in the presence of convective flows. We go for stochastic features in our model and versions of operations to have more accurate predictions and to know more new features in the behaviour of these dynamic systems.

In this work, we shall discuss the theoretical foundation of the fluid dynamics of Casson fluids, the effect of varying thermal conductivity and mass diffusivity in the problems, and introduce some stochastic numerical algorithms for solving such complex flow systems.

The primary contribution of this work is the suggestion of a stochastic numerical scheme for the solution of stochastic partial differential equations. Another method exists in the literature for solving stochastic partial differential equations. That scheme is called the Maruyama method and can be used to solve stochastic equations that appear in fluid dynamics with the variation of Wiener process terms. The Euler-Maruyama method extends the more common forward Euler approach for stochastic differential equations. The Matlab commands generate random numbers from a Normal distribution with a mean of zero and a standard deviation that determines the time step size for the Wiener process term in the scheme.

This contribution’s stochastic numerical approach can be utilized to solve partial differential equations. The scheme is based on two steps. A partial differential equation’s solution can be predicted in the first step, the predictor stage. The second stage is the corrector stage, which finds the partial differential equation’s solution. But these two stages only find the solution for the deterministic model. The scheme for finding the solution of the stochastic differential equation will be proposed later. For proposing a scheme for a deterministic equation, consider the deterministic equation as follows: ∂ v ∂ t = G v , ∂ v ∂ x , ∂ v ∂ y , ∂ 2 v ∂ y 2 (1)

Let the first stage of the scheme be expressed as: v ¯ i , j n + 1 = v i , j n + Δ t   ∂ v ∂ t i , j n (2)

Where Δ t is the time step size.

The second stage of the scheme is expressed as: v i , j n + 1 = 1 5 4 v i , j n + v ¯ i , j n + 1 + Δ t a ∂ v ∂ t i , j n + b ∂ v ¯ ∂ t i , j n + 1 (3) where a and b will be determined later.

Now, substitute Eq. 2 into Eq. 3, which yields. v i , j n + 1 = 1 5 4 v i , j n + v i , j n + Δ t ∂ v ∂ t i , j n + Δ t a   ∂ v ∂ t i , j n + b   ∂ v ∂ t i , j n + b Δ t   ∂ 2 v ∂ t 2 i , j n (4)

Expanding v i , j n + 1 using Taylor series expansion v i , j n + 1 = v i , j n + Δ t   ∂ v ∂ t i , j n + Δ t 2 2 ∂ 2 v ∂ t 2 i , j n + O Δ t 3 (5)

Substituting Eq. 5 into Eq. 4 gives v i , j n + Δ t ∂ v ∂ t i , j n + Δ t 2 2 ∂ 2 v ∂ t 2 i , j n = v i , j n + 1 5 Δ t ∂ v ∂ t i , j n + Δ t a   ∂ v ∂ t i , j n + b   ∂ v ∂ t i , j n + b Δ t   ∂ 2 v ∂ t 2 i , j n (6)

Evaluating the coefficients of ∂ v ∂ t i , j n and ∂ 2 v ∂ t 2 i , j n on both sides of Eq. 6 that yields 1 = 1 5 + a + b 1 2 = b (7)

Solving Equation 7 gives a = 3 10 , b = 1 2 (8)

Therefore, the time discretization of Eq. 1 is v ¯ i , j n + 1 = v i , j n + Δ t G v i , j n , ∂ v ∂ x i , j n , ∂ v ∂ y i , j n , ∂ 2 v ∂ y 2 i , j n (9) v i , j n + 1 = 1 5 ( 4 v i , j n + v ¯ i , j n + 1 ) + Δ t a   G v i , j n , ∂ v ∂ x i , j n , ∂ v ∂ y i , j n , ∂ 2 v ∂ y 2 i , j n + b   G v ¯ i , j n + 1 , ∂ v ¯ ∂ x i , j n + 1 , ∂ v ¯ ∂ y i , j n + 1 , ∂ 2 v ¯ ∂ y 2 i , j n + 1 (10)

Now consider the partial differential equations as d v = G v , ∂ v ∂ x , ∂ v ∂ y , ∂ 2 v ∂ y 2 d t + σ f v d W t (11)

Its I ^ t o integral form is given as v n + 1 = v n + ∫ t n t n + 1 G d t + ∫ t n t n + 1 σ f v d W t (12)

Using Taylor series expansion for f v as f v = f v n + Δ t f ′ v n + Δ t 2 2 f ″ v n + O Δ t 3 (13)

So, the last term in Eq. 12 can be expressed as ∫ t n t n + 1 f v d W = ∫ t n t n + 1 f v n + Δ t f ′ v n + Δ t 2 2 f ″ v n d W = f v = f v n Δ W + Δ t f ′ v n Δ W + Δ t 2 2 f ″ v n Δ W (14)

Therefore, the proposed stochastic numerical scheme for time discretization Eq. 11 is v ¯ i , j n + 1 = v i , j n + Δ t G v i , j n , ∂ v ∂ x i , j n , ∂ v ∂ y i , j n , ∂ 2 v ∂ y 2 i , j n (15) v i , j n + 1 = 1 5 4 v i , j n + v ¯ i , j n + 1 + Δ t { 3 10 G v i , j n , ∂ v ∂ x i , j n , ∂ v ∂ y i , j n , ∂ 2 v ∂ y 2 i , j n + 1 2 G ( v ¯ i , j n + 1 , ∂ v ¯ ∂ x i , j n + 1 , ∂ v ¯ ∂ y i , j n + 1 , ∂ 2 v ¯ ∂ y 2 i , j n + 1 ) } + σ f v n Δ W + σ Δ t f ′ v n Δ W + σ 2 Δ t 2 f ″ v n Δ W (16) where Δ W is approximated as a normal distribution with mean 0 and standard deviation Δ t i . e . , Δ W ∼ N 0 , Δ t

Let f v = v and G = β 1 ∂ q v ∂ x + β 2 ∂ q v ∂ y + β 3 ∂ q 2 v ∂ y 2 then the fully discretized scheme is given as v ¯ i , j n + 1 = v i , j n + Δ t β 1 δ x v i , j n + β 2 δ y   v i , j n + β 3 δ y 2 v i , j n q + 1 2 (17) v i , j n + 1 = 1 5 4 v i , j n + v ¯ i , j n + 1 + Δ t { 3 10 β 1 δ x v i , j n + β 2 δ y v i , j n + β 3 δ y 2 v i , j n q + 1 2 + 1 2 ( β 1 δ x v ¯ i , j n + 1 + β 2 δ y   v ¯ i , j n + 1 + β 3 δ y 2 v ¯ i , j n + 1 1 + q 2 ) + σ v i , j n Δ W + σ Δ t Δ W } (18) where δ x v i , j n = v i + 1 , j n − v i − 1 , j n 2 Δ x , δ y v i , j n = v i , j + 1 n − v i , j − 1 n 2 Δ y , δ y 2 v i , j n = v i , j + 1 n − 2 v i , j n + v i , j − 1 n Δ y 2

Theorem 1

The proposed numerical schemes (17) and (18) are consistent in the mean square sense.

Proof. Let P be a smooth function, then. L P i n = P n + 1 Δ t , i Δ x , j Δ y − P n Δ t , i Δ x , j Δ y − β 1 ∫ n Δ t n + 1 Δ t P x s , i Δ x , j Δ y d S − β 2 ∫ n Δ t n + 1 Δ t P y s , i Δ x , j Δ y d S − β 3 ∫ n Δ t n + 1 Δ t P y y s , i Δ x , j Δ y d S − σ ∫ n Δ t n + 1 Δ t P s , i Δ x , j Δ y ) d W ( S (19)

Now, combining both stages of the schemes gives the following operator L i n P = P n + 1 Δ t , i Δ x , j Δ y − P n Δ t , i Δ x , j Δ y − Δ t [ β 1 4 Δ x P n Δ t , i + 1 Δ x , j Δ y − P n Δ t , i − 1 Δ x , j Δ y + β 2 4 Δ y P n Δ t , i Δ x , j + 1 Δ y − P n Δ t , i Δ x , j − 1 Δ y + β 3 2 Δ y 2 ( P n Δ t , i Δ x , j + 1 Δ y − 2 P n Δ t , i Δ x , j Δ y + P n Δ t , i Δ x , j − 1 Δ y ) ] − Δ t [ β 1 4 Δ x ( P ¯ n + 1 Δ t , i + 1 Δ x , j Δ y − P ¯ n + 1 Δ t , i − 1 Δ x , j Δ y ) + β 2 4 Δ y ( P ¯ n + 1 Δ t , i Δ x , j + 1 Δ y − P ¯ n + 1 Δ t , i Δ x , j − 1 Δ y ) + β 3 2 Δ y 2 ( P ¯ n + 1 Δ t , i Δ x , j + 1 Δ y − 2 P ¯ n + 1 Δ t , i Δ x , j Δ y + P ¯ n + 1 Δ t , i Δ x , j − 1 Δ y ) ] − σ P n Δ t , i Δ x , j Δ y W n + 1 Δ t − W n Δ t − σ Δ t W n + 1 Δ t − W n Δ t (20) where P ¯ n + 1 Δ t , i Δ x , j Δ y = P n Δ t , i Δ x , j Δ y + Δ t { β 1 2 Δ x P n Δ t , i + 1 Δ x , j Δ y − P n Δ t , i − 1 Δ x , j Δ y + β 2 2 Δ y P n Δ t , i Δ x , j + 1 Δ y − P n Δ t , i Δ x , j − 1 Δ y + β 3 Δ y 2 P n Δ t , i Δ x , j + 1 Δ y − 2 P n Δ t , i Δ x , j Δ y + P n Δ t , i Δ x , j − 1 Δ y }

The mean square of the scheme is written as: E L P i n − L i n P 2 = E | − β 1 ∫ n Δ t n + 1 Δ t P x s , i Δ x , j Δ y d S − β 2 ∫ n Δ t n + 1 Δ t P y s , i Δ x , j Δ y d S − β 3 ∫ n Δ t n + 1 Δ t P y y s , i Δ x , j Δ y d S − σ ∫ n Δ t n + 1 Δ t P s , i Δ x , j Δ y d W S + β 1 Δ t 4 Δ x P n Δ t , i + 1 Δ x , j Δ y − P n Δ t , i − 1 Δ x , j Δ y + β 2 Δ t 4 Δ y P n Δ t , i Δ x , j + 1 Δ y − P n Δ t , i Δ x , j − 1 Δ y + β 3 Δ t 2 Δ y 2 ( P n Δ t , i Δ x , j + 1 Δ y − 2 P n Δ t , i Δ x , j Δ y + P n Δ t , i Δ x , j − 1 Δ y ) + β 1 Δ t 4 Δ x P ¯ n + 1 Δ t , i + 1 Δ x , j Δ y − P ¯ n + 1 Δ t , i − 1 Δ x , j Δ y + β 2 Δ t 4 Δ y P ¯ n + 1 Δ t , i Δ x , j + 1 Δ y − P ¯ n + 1 Δ t , i Δ x , j − 1 Δ y + β 3 Δ t 2 Δ y 2 ( P ¯ n + 1 Δ t , i Δ x , j + 1 Δ y − 2 P ¯ n + 1 Δ t , i Δ x , j Δ y + P ¯ n + 1 Δ t , i Δ x , j − 1 Δ y ) + σ Δ t + P n Δ t , i Δ x , j Δ y × W n + 1 Δ t − W n Δ t | 2 (21)

Equation 21 can be written as: E L P i n − L i n P 2 ≤ 2 β 1 2 E | ∫ n Δ t n + 1 Δ t − P x s , i Δ x , j Δ y d S + Δ t 4 Δ x P n Δ t , i + 1 Δ x , j Δ y − P n Δ t , i − 1 Δ x , j Δ y + Δ t 4 Δ x P ¯ n + 1 Δ t , i + 1 Δ x , j Δ y − P ¯ n + 1 Δ t , i − 1 Δ x , j Δ y | 2 + 2 β 2 2 E ∫ n Δ t n + 1 Δ t − P y s , i Δ x , j Δ y d S + Δ t 4 Δ y P n Δ t , i Δ x , j + 1 Δ y − P n Δ t , i Δ x , j − 1 Δ y + Δ t 4 Δ y P ¯ n + 1 Δ t , i Δ x , j + 1 Δ y − P ¯ n + 1 Δ t , i Δ x , j − 1 Δ y | 2 + 2 β 3 2 E | ∫ n Δ t n + 1 Δ t − P y y s , i Δ x , j Δ y d S + Δ t 2 Δ y 2 P n Δ t , i Δ x , j + 1 Δ y − 2 P n Δ t , i Δ x , j Δ y + P n Δ t , i Δ x , j − 1 Δ y + Δ t 2 Δ y 2 P ¯ n + 1 Δ t , i Δ x , j + 1 Δ y − 2 P ¯ n + 1 Δ t , i Δ x , j Δ y + P ¯ n + 1 Δ t , i Δ x , j − 1 Δ y | 2 + 2 σ 2 E | ∫ n Δ t n + 1 Δ t − P s , i Δ x , j Δ y d W S Δ t + Q n Δ t , i Δ x , j Δ y W n + 1 Δ t − W n Δ t | 2 (22)

Now, utilizing the result E ∫ n Δ t n + 1 Δ t − P s , i Δ x , j Δ y − Δ t + Q n Δ t , i Δ x , j Δ y d W ( S 2 ≤ Δ t ∫ n Δ t n + 1 Δ t E − P s , i Δ x , j Δ y − Δ t + Q n Δ t , i Δ x , j Δ y 2 d S (23)

Thus by applying l i m i t as Δ x → 0 , Δ y → 0 , Δ t → 0 and n Δ t , i Δ x , j Δ y → t , x , y , the mean square error approaches zero. i . e . E L P i n − L i n P 2 → 0 (24)

So, the proposed scheme is consistent.

Theorem 2

The proposed numerical scheme is conditionally stable.

Proof: The stability of the proposed scheme will be analyzed using Fourier series analysis and mean square sense. The Fourier series analysis for the classical model will be applied, and then stability conditions in the mean square sense will be employed. The Fourier series analysis requires some transformations when finding stability conditions of finite difference schemes. The transformation reduces the difference equation into trigonometric equations, and the stability condition will be determined later. For applying a Taylor series analysis for scheme (17) and (18), the following transformations will be applied v ¯ i , j n + 1 = Q ¯ n + 1 e i I ψ 1 e j I ψ 2 , v i , j n = Q n e i I ψ 1 e j I ψ 2 v i ± 1 , j n = Q n e i ± 1 I ψ 1 e j I ψ 2 , v i , j ± 1 n = Q n e i I ψ 1 e j ± 1 I ψ 2 v ¯ i ± 1 , j n + 1 = Q ¯ n + 1 e i ± 1 I ψ 1 e j I ψ 2 , v ¯ i , j ± 1 n + 1 = Q ¯ n + 1 e i I ψ 1 e j ± 1 I ψ 2 (25)

Applying some of the transformations from Eq. 25 in the first stage of scheme (17) yields. Q ¯ n + 1 e i I ψ 1 e j I ψ 2 = Q n e i I ψ 1 e j I ψ 2 + Δ t β 1 e i + 1 I ψ 1 e j I ψ 2 − e i − 1 I ψ 1 e j I ψ 2 2 Δ x Q n + β 2 e i I ψ 1 e j + 1 I ψ 2 − e i I ψ 1 e j − 1 I ψ 2 2 Δ y Q n + β 3 e i I ψ 1 e j + 1 I ψ 2 − 2 e i I ψ 1 e j I ψ 2 + e i I ψ 1 e j − 1 I ψ 2 Δ y 2 Q n (26)

Upon dividing both sides of Eq. 27 by e i I ψ 1 e j I ψ 2 , it yields Q ¯ n + 1 = Q n + Δ t β 1 2 Δ x e I ψ 1 − e − I ψ 1 + β 2 2 Δ y e I ψ 2 − e − I ψ 2 + β 3 Δ y 2 e I ψ 2 − 2 + e − I ψ 2 Q n (27)

Using trigonometric identities in Eq. 27 it yields Q ¯ n + 1 = Q n + Δ t β 1 Δ x I s i n ψ 1 + β 2 Δ y I s i n ψ 2 + 2 β 3 Δ y 2 cos ⁡ ψ 2 − 1 Q n (28)

Re-write Eq. 28 as: Q ¯ n + 1 = Q n + c 1 I s i n ψ 1 + c 2 I s i n ψ 2 + c 3 cos ⁡ ψ 2 − 1 Q n (29) where c 1 = β 1 Δ t Δ x , c 2 = β 2 Δ t Δ y , c 3 = 2 β 3 Δ t Δ y 2

Similarly, employing some of the transformation from Eq. 25 into the second stage of the scheme and ignoring the non-homogeneous part in Eq. 18 gives Q n + 1 = 1 5 Q n + Q ¯ n + 1 + 3 10 c 1 I s i n ψ 1 + c 2 I s i n ψ 2 + c 3 cos ⁡ ψ 2 − 1 Q n + 1 2 c 1 I s i n ψ 1 + c 2 I s i n ψ 2 + c 3 cos ⁡ ψ 2 − 1 Q ¯ n + 1 + σ Q n Δ W (30)

Substituting Eq. 29 into Eq. 30 yields Q n + 1 = Q n + 1 2 c 1 I s i n ψ 1 + c 2 I s i n ψ 2 + c 3 cos ⁡ ψ 2 − 1 Q n + 1 2 c 1 I s i n ψ 1 + c 2 I s i n ψ 2 + c 3 cos ⁡ ψ 2 − 1 c 1 I s i n ψ 1 + c 2 I s i n ψ 2 + c 3 cos ⁡ ψ 2 − 1 Q n + σ Δ W Q n (31)

Re-write Eq. 31 as Q n + 1 = Q n + 1 2 Z Q n + 1 2 Z 1 + Z Q n + σ Δ W Q n (32) where Z = c 1 I s i n ψ 1 + c 2 I s i n ψ 2 + c 3 cos ⁡ ψ 2 − 1

Equation 32 can be re-written as Q n + 1 = a ¯ + I b ¯ Q n + σ Δ W Q n (33)

The amplification factor is written as Q n + 1 Q n 2 = a ¯ + σ Δ W 2 + b ¯ 2 (34) where a ¯ = R e   Z + 1 2 R e   Z 2 − I m   Z 2 and b ¯ = I m   Z + R e   Z   I m   Z

Applying expected value on both sides of Eq. 33 yields E Q n + 1 Q n 2 = E a ¯ 2 + b ¯ 2 + 2 σ a ¯ E Δ W + σ 2 E Δ W 2 (35)

Since E Δ W = 0 , a n d   E Δ W 2 = Δ t

Therefore, Eq. 35 yields E Q n + 1 Q n 2 = a ¯ 2 + b ¯ 22 + σ 2 Δ t (36)

Now if a ¯ 2 + b ¯ 22 ≤ 1 and let λ = σ 2 then Eq. 37 can be re-written as E Q n + 1 Q n 2 = 1 + λ Δ t (37)

Thus, the proposed stochastic numerical scheme with non-homogenous parts is conditionally stable in the mean square sense.

Below, we present a comprehensive analysis of the advantages and disadvantages of the proposed scheme.

Consider the non-Newtonian, unsteady, laminar, and incompressible fluid flow over the sheet. The plate’s abrupt motion induces fluid flow toward the positive x * -axis, where the x * -axis represents the horizontal direction, and the y * -axis is perpendicular to it. The stretching velocity of the plate is represented by u w . A uniform electric field E ⇀ = 0 , 0 , − E ∘ and transverse magnetic field B ⇀ = 0 , B ∘ , 0 are applied, and the fluid is electrically conducting. The electric and magnetic fields follow Ohm’s rule, but the electric field is stronger. For the moment, disregard the Hall effect and the induced magnetic field. Chemical reactions, frictional heating, and viscous dissipation are some of the flow characteristics taken into account. Under the assumption of boundary theory over a flat plate, the governing equations are expressed as: ∂ u * ∂ x * + ∂ v * ∂ y * = 0 (38) ∂ u * ∂ t * + u * ∂ u * ∂ x * + v * ∂ u * ∂ y * = ν 1 + 1 β ∂ 2 u * ∂ y * 2 + σ ρ E ∘ B ∘ − B ∘ 2 u * + g β T T − T ∞ + β C C − C ∞ (39) ∂ T ∂ t * + u * ∂ T ∂ x * + v * ∂ T ∂ y * = 1 ρ C p ∂ ∂ y * k T ∂ T ∂ y *   + ν C p 1 + 1 β ∂ u * ∂ y * 2 + σ ρ C p u B ∘ − E ∘ 2 (40) ∂ C ∂ t * + u * ∂ C ∂ x * + v * ∂ C ∂ y * = ∂ ∂ y * D C ∂ C ∂ y * − k r C − C ∞ (41)

With the following initial and boundary conditions u * = 0 , v * = 0 , T = 0 , C = 0 , w h e n   t * = 0 u * = u w , v * = 0 , T = T w , C = C w , w h e n   y * = 0 u * → 0 , T → T ∞ , C → C ∞ , w h e n   y * → ∞ u * = 0 , v * = 0 , T = 0 , C = 0   w h e n   x * = 0 (42) where k T = k ∞ 1 + ε θ and D C = D ∞ 1 + ε ¯ ϕ and σ is electrical conductivity, ρ is the density of the fluid, C p is specific heat capacity, β is the Casson parameter, g is the gravity, B T is the coefficient of thermal expansion and β C is the coefficient of solutal expansion and k r is reaction rate. The transformations u = u * u w , v = v * u w , t = t * u w L , x = x * L , y = y * L , θ = T − T ∞ T w − T ∞ , ϕ = C − C ∞ C w − C ∞ (43)

When applied to Eqs. 38–42 reduces them to following dimensionless equations ∂ u ∂ x + ∂ v ∂ y = 0 (44) ∂ u ∂ t + u ∂ u ∂ x + v ∂ u ∂ y = 1 R e 1 + 1 β ∂ 2 u ∂ y 2 + H a 2 R e E 1 − u + G γ T R e 2 θ + G γ C R e 2 ϕ (45) ∂ θ ∂ t + u ∂ θ ∂ x + v ∂ θ ∂ y = ε P r 1 R e ∂ θ ∂ y 2 + 1 P r 1 R e 1 + ε θ ∂ 2 u ∂ y 2 + E C H o 2 R e u − E 1 2 + E C R e 1 + 1 β ∂ u ∂ y 2 (46) ∂ ϕ ∂ t + u ∂ ϕ ∂ x + v ∂ ϕ ∂ y = ε ¯ S c R c ∂ ϕ ∂ y 2 + 1 S c 1 R e 1 + ε ¯ ϕ ∂ 2 ϕ ∂ y 2 − γ ϕ (47)

Subject to the dimensionless boundary and initial conditions u = 0 , v = 0 , θ = 0 , ϕ = 0   w h e n   t = 0 u = 1 , v = 0 , θ = 1 , ϕ = 1   w h e n   y = 0 u → 0 , θ → 0 , ϕ → 0   w h e n   y → ∞ u = 0 , v = 0 , θ = 0 , ϕ = 1   w h e n   x = 0 (48) where H a is Hartmann’s number, E 1 is used for local electric parameters, G γ T is thermal Grashof number, G γ C solutal Grashof number, E C Eckert number, R e is Reynolds number, P r is Prandtl number, S c is Schmidt number and γ dimensionless reaction rate parameter, and these are defined as H a = σ ρ ν B ∘ L , E 1 = E ∘ B ∘ u w , G γ T = L 3 g β T T w − T ∞ ν 2 , G Y C = L 3 g β C C w − C ∞ ν 2 , E C = u w 2 C p T w − T ∞ , R e = L u w ν , P r = ρ C p ν k ∞ , S c = D ∞ ν , γ = L k γ u w