Advancements in Numerical Solutions: Fractal Runge-Kutta Approach to Model Time-Dependent MHD Newtonian Fluid with Rescaled Viscosity on Riga Plate

Nomenclature

1  Introduction

Fluid dynamics have long enthralled scientific investigation due to its complicated and sophisticated interplay of forces. Researchers have recently focused more on solving the problems of time-dependent phenomena in fluid dynamics, thanks to improvements in numerical approaches and a spike in processing power [1]. Fractal time-dependent problems are incredibly fascinating since they present unique difficulties for numerical analysis due to their complexity and self-repeating patterns. A unique computational approach is presented in this study to tackle the difficulties of fractal time-dependent issues and contribute to the improvement of numerical solutions in fluid dynamics [2]. We are particularly interested in the area where Newtonian fluid dynamics meets magnetohydrodynamics (MHD), where the behaviour of time-dependent fluids affected by magnetic fields is crucial. Another new development in geometry is ‘fractal geometry’, which helps to create a masterpiece of the flow pattern and its general attributes on several scales [3]. This can significantly improve the results and the way of modelling behind employing complexities and methods of fractal view.

The purpose of the current study was to verify the complexity of flow by comparing it to other fractals. Fractals in mathematics are geometric appearances that can be divided into pieces and cut out, where each piece’s reduction is that very essence known as fractality. We need to adjust these linear and non-linear differential equations that give solutions of flow fields with a set of calculations so it can quickly figure out the different patterns of flow fields over several times and space scales. We use the Fractal Runge-Kutta Method to solve the problem. It is an efficient way to find an approximate answer, whereas when we use the method in flow, it successfully displays the flow details at different times and scales. When we try to study the dynamics of a time-dependent fractal MHD Newtonian fluid with rescaled viscosity over the Riga plate, a moving surface is used as the wall [4]. Since identifying the flow detail and heat transfer is not a simple problem, we also needed the terminal of mathematical methods to derive the results. Therefore, we used a fractal method to model the problem. The comparison with fractals as and when needed can help us determine the flow pattern’s complexity. In the modern techniques of geometry or the field of fluid flow analysis, the employment of fractal theory can improve the predictability of the pattern of flow behaviour. In the study, we see the use of geometry in the case of a dynamics problem [5].

Flows similar at different time scales can reflect the complexity of the fractal structure of a flow. There are various ways to make a computer do fractal analysis, and various computer programs are available on the internet to help us perform that task on our own. One such famous fractal is far from us to find a set of complex roots of a polynomial. We needed an efficient way to get a target result with a close approach; thus, we used a completely modified approach, Runge-Kutta. We can infinitely approximate the error of the flow field into a particular area and then fix those results and validate them from our fractal comparison.

The Fractal Runge-Kutta Method is a version of the Runge-Kutta scheme and is a very efficient way to solve differential equations that show fractal characteristics. Traditionally, the equations to be solved are essentially third-, fourth-, and fifth-order differential equations. Hence, a suitable fractal point is to be achieved, and then the differential equations are solved to get a particular result. That’s where Fractal Runge Kutta differs from traditional. It comes up with a target point and the error value to get an exact result value that is as close as possible to the append error. And that’s what we needed to solve the problem. It allows the accurate portrayal of physical phenomena. To summarise, Fractals can assist us in many fractals when we get stuck. Solving difficult problems yields constants-based fractal calculations that allow us to draw quick conclusions and retain seemingly impossible answers. These fractals would account for many mathematical interpretations while we analyze the problem through different patterns. When we try to see what is happening in the flow (when the flow field turns out to be very subtle), we find connecting these Roman coefficients to form a Pythagorean Theorem and then adjusting the results as much as possible.

A boundary layer in fluid dynamics is a thin layer that has been researched in theory in fluid dynamics. It’s a thin layer of each phase or a boundary layer in each phase. It has been around for more than a century. It started when Ludwig Prandtl developed it in 1904. The boundary layer has powerful effects on the motion of fluids and represents their characteristics. Hence, studying a Riga plate is an essential dynamic phenomenon [6]. It can be stated as a fluid layer that flows directly over a surface, including shear flow. Researchers have conducted studies to explore the nature of the flow of the layers in various forms of geometry, elongated shape, and horizontal plates [7], Read Less vertical plates, being infinite in length [8], stretching sheets [9], and Read Less wedges [10]. A novel device called the Riga plate has been developed to minimize skin friction in movement [11]. The electromagnetic properties of the Riga plate can be utilized to study the effects of thermal radiation and viscous dissipation when the plate is submerged in a setup by positioning electrodes on it. Research on this geometry has lately gathered steam in the 21st Century, despite Gailitis and Lielausis introducing this plate in 1961 [12]. Using the Cattaneo Christov model [13], we studied the impact of radiation and chemical reactions on Darcy Forchheimer flow, focusing on a Williamson fluid passing over a Riga plate. Later, Reyaz et al. [14] used both numerical methods to analyze the data.

In a study [14], the Zakians method was employed to analyze the flow behaviour near a Riga plate under the influence of a magnetic field. The findings suggested that the presence of the Riga plate impacts the velocity. Particle deposition effects in Casson-type hybrid nanofluid flow over a Riga plate were investigated in another study by [15] utilizing techniques like the Runge-Kutta Fehlberg method. Additionally, using analysis [16]. They have examined the stagnation point flow of a hybrid nanofluid across a Riga plate. Furthermore, in their paper, Loganathan et al. [17] examined how variable chemical reactions affect Casson fluid flow over a Riga plate in a medium. The buoyancy-driven Casson fluid flow from a vertical wavy surface in the presence of a magnetic field is studied numerically in the proposed study [18], which also provides a quantitative and qualitative assessment of the underlying heat transport in the free convective regime.

Computer networking, dynamics, control theory, and biology have all seen a rise in the use of fractional calculus. Analyzing memory effects is made easier using this tool. Fluid problems are now being approached using methods, a new trend [19]. Different mathematicians have defined derivatives with and without kernels [20]. Atangana Baleanu derivative is one example of such a derivative [21]. When dealing with nonlocal behaviour and equations with non-singular kernels, this derivative proves to be helpful [22]. It assists in addressing real-world issues such as heat and mass transfer, electrical circuits, image processing, and the adverse effects of chemotherapy in cancer treatment [23]. The Casson fluid flow alongside a plate was studied by Tassaddiq et al. [24] using a fractional approach. Using the Laplace transform and Zakians algorithm, they solved transformed equations.

Shah and colleagues examined the Atangana-Baleanu derivative’s utilization in the flow analysis between plates [25] and cylinders [26]. The study in [27] investigated the flow behaviour of nanofluids containing spherical-shaped Molybdenum Disulphide nanoparticles using an AB derivative. Saqib et al. [28] examined the CNT nanofluid’s channel flow properties with CMC as the base fluid employing an AB derivative. This numerical study [29] investigated the fractional viscoelastic fluid model with unstable convection on an inclined plane through a Forchheimer medium. This paper [30] focuses on developing a fractional explicit-implicit numerical method for solving time-dependent partial differential equations. In their work, Nawaz et al. [31] proposed a third-order numerical approach for solving fractional partial differential equations. The initial explicit stage can rapidly converge, while the subsequent implicit stage is responsible for expanding the stability region.

Numerous research efforts have been devoted to exploring the practical applications of nanofluids in heat transfer studies. Nanofluids are essentially solutions comprised of a fluid infused with nanoparticles synthesized in laboratories to create an innovative solution that enhances overall fluid behaviour for improved outcomes. While work on these fluids has been ongoing since the 1990s, Choi first introduced and coined the term “nanofluids” in 1995 [32]. Using far nanofluids has yielded promising outcomes in various disciplines, from engineering to healthcare. The consequences of nanofluid flow in geometries have been the subject of much mathematical discussion; these findings have significant practical ramifications.

Offer insights into the utilization of nanoparticles in the base fluid. Magnetite nanoparticles exhibit potential in applications and research [33]. They prove helpful for diagnosing and treating diseases. Due to their potential in applications, researchers have devoted significant attention to magnetite nanoparticles in their studies for several decades. The features of these nanoparticles make them widely used in treatment. The effect of conductivity on magnetite nanofluid-based sheets subjected to shrinkage and stretching was studied by [34]. In their paper, Seth et al. [35] discussed ferrofluid flow over a plate with oscillating behaviour.

In this context, while traditional Runge-Kutta (RK) methods effectively solve equations, they often face challenges when capturing the complex and non-linear dynamics of time fractal MHD Newtonian fluid with rescaled viscosity flows over Riga plates. This is where Fractal Runge-Kutta (FRK) methods emerge as an alternative. FRK methods utilize the concept of calculus, extending the idea of differentiation and integration beyond numbers. Complex physical systems frequently exhibit fractal time’s complexity and irregular structure, which exemplifies it. By incorporating these derivatives into the RK framework, FRK methods demonstrate improved accuracy and stability, especially when dealing with highly non-linear and singular problems.

The previous studies on flow and heat transfer on stretching sheets have captured the attention of scholars across fields. It holds significance in polyamide production due to advancements in polymers. With the idea of moving planes and velocities that change the distance from fixed locations on a sheet in a linear fashion, Crane [36] investigated these concepts. Many researchers have recently employed Newtonian fluids for heat and mass transfer [37,38]. Nadeem et al. [39] investigated second-grade fluids’ time-dependent stretching behaviour. The paper by Majeed et al. [40] suggested using the suction of stretching surfaces to study the flow of non-Newtonian fluids. The References [41,42] described a new class of exact solutions to the Oberbeck-Boussinesq equations that describe an incompressible fluid.

1.1 Novelty of the Study

This study is unique because it takes an interdisciplinary approach, applies the Fractal Runge-Kutta Method, develops an efficient computational scheme, investigates fractal time-dependent fluid dynamics, and uncovers valuable insights into fluid dynamics phenomena. All of these things add up to a more advanced numerical solution to complicated fluid dynamic problems, which in turn opens up new possibilities for study and development in the area:

1. A new computational approach is used to solve fractal time-dependent fluid dynamics problems. The Fractal Runge-Kutta Method models time-dependent MHD Newtonian fluid with rescaled viscosity on Riga plates.

2. A three-stage computational approach using three-time levels is proposed. This comprehensive approach captures the complexity of fractal time-dependent problems to provide a sophisticated knowledge of fluid dynamics.

3. Fourier series analysis of scalar issues proves the computing scheme’s stability. Moreover, a convergence analysis for a time-dependent fractal partial differential equations system ensures the numerical method’s dependability and correctness.

4. The computational approach is applied to a dimensionless fractal model of incompressible, unstable, laminar Newtonian fluid with rescaled viscosity over flat and oscillatory Riga plates. For practicality, space- and temperature-dependent heat sources are included in the numerical simulations.

5. We present a powerful computational technique for fractal time-dependent fluid dynamics problems. We ensure numerical solution accuracy and efficiency using a scheme with three-time levels and rigorous stability analysis.

1.2 Potential Applications

The Riga plate geometry was selected for our investigation due to theoretical and practical considerations of fluid dynamics. What follows is an outline of our decision-making procedure, along with some examples of its potential applications:

1. Theoretical Considerations: The Riga plate design is straightforward and analytically understandable for fluid flow studies. Its primary and adjustable geometry allows for studying fundamental fluid dynamics topics, making it a popular choice for theoretical and computational fluid mechanics investigations. We can use the Riga plate’s simple geometry to learn about time-varying fluid dynamics in controlled conditions.

2. Practical Applications: Despite its theoretical simplicity, the Riga plate geometry has extensive practical applications across various industries. For example:

a. Heat Exchangers: The Riga plate geometry finds utility in heat exchangers, which transfer heat from one fluid to another using a series of flat or oscillating plates. A comprehensive comprehension of the fluid dynamics across Riga plates is necessary for systems to achieve optimal heat transfer efficiency.

b. Aerospace Engineering: Regarding aerospace applications, the Riga plate geometry is a common feature in the design of airfoils and wings. The fluid flow characteristics over these surfaces are essential in determining the aerodynamic performance. Studying fluid dynamics over Riga plates can provide valuable insights into airflow behaviour, which can help improve efficiency and stability in aerospace systems.

Therefore, the article focuses on advancing and implementing a particular FRK technique to investigate the time fractal MHD Newtonian fluid dispersion with rescaled viscosity flows on Riga plates. We expect that this method will offer valuable information regarding the correlation between the Newtonian characteristics of the fluid, magnetic fields, and fractal time scales. These insights can be applied in various sectors, including optimizing heat transfer processes, regulating flow instabilities, and creating innovative material coatings.

2  The Fractal Numerical Scheme

The fractal scheme is a numerical method consisting of three stages that can be used to discretize time-dependent PDEs. Any finite difference discretization or other scheme can be employed for spatial discretization. Before proposing a scheme, the whole domain is divided into small subintervals with endpoints called grid points. The solution will be found at each grid point using the proposed scheme except at those where boundary conditions are given. The main advantage of using the proposed scheme is that it does not require linearization or any other iterative scheme. The scheme is employed to follow time fractal partial differential equations.

∂v∂tα=G(v,∂v∂x,∂v∂y,∂2v∂y2)(1)

where 0<α≤1. Subject to initial and boundary conditions given as:

v(0,x,y)=0v(t,0,y)= α1  v(t,L,y)= α2,v(t,x,0)=α3v(t,x,L)=α4}(2)

where αi, i=1,2,3,4 are constant.

The first stage of the scheme is:

v¯i,jn+1=vi,jn+Δt1∂v∂tα|i,jn(3)

The first stage (1) finds the solution of Eq. (1) at an arbitrary time level. This stage (1) can also be called as predictor stage. The second stage of the proposed scheme can be written as:

v=i,jn+1=12(vi,jn+v¯i,jn+1)+Δt2∂v¯∂tα|i,jn+1(4)

This stage (4) is also a predictor stage that finds the solution of the fractal partial differential equation at an arbitrary time level. It also utilizes the solution computed from the first predictor stage. The third stage of the proposed scheme can be articulated as:

vi,jn+1=avi,jn+bv¯i,jn+cv=i,jn+1+dΔt∂v=∂tα|i,jn+1(5)

where a,b,c, and d are parameters to be determined.

Now, substituting the first predictor stage (3) into the second predictor stage (4) yields:

v=i,jn+1=vi,jn+12Δt1∂v∂tα|i,jn+Δt2∂v∂tα|i,jn+Δt1Δt2(∂∂tα)2vi,jn(6)

By using Eqs. (3) and (6) into corrector stage (5), it yields:

vi,jn+1=avi,jn+bvi,jn+bΔt1∂v∂tα|i,jn+cvi,jn+c2Δt1∂v∂tα|i,jn+cΔt2∂v∂tα|i,jn+cΔt1Δt2(∂∂tα)2vi,jn+dΔt{∂v∂tα|i,jn+12Δt1(∂∂tα)2vi,jn+Δt2(∂∂tα)2vi,jn+Δt1Δt2(∂∂tα)3vi,jn}(7)

Expanding vi,jn+1 using fractal Taylor series yields:

vi,jn+1=vi,jn+Δt∂v∂tα|i,jn+(Δt)22(∂∂tα)2vi,jn+(Δt)36(∂∂tα)3vi,jn+O((Δt)4)(8)

By substituting fractal Taylor series expansion of vi,jn+1 (8) into Eq. (7), it is obtained:

vi,jn+Δt∂v∂tα|i,jn+(Δt)22(∂∂tα)2vi,jn+(Δt)36(∂∂tα)3vi,jn=avi,jn+bvi,jn+bΔt1∂v∂tα|i,jn+cvi,jn+c2Δt1∂v∂tα|i,jn+cΔt2∂v∂tα|i,jn+cΔt1Δt2(∂∂tα)2vi,jn+dΔt{∂v∂tα|i,jn+12Δt1(∂∂tα)2vi,jn+Δt2(∂∂tα)2vi,jn+Δt1Δt2(∂∂tα)3vi,jn}(9)

By equating coefficients of vi,jn,∂v∂tα|i,jn,(∂∂tα)2vi,jn and (∂∂tα)3vi,jn on both sides of Eq. (9), yields:

a+b+c=1bΔt1+c2Δt1+cΔt2+dΔt=ΔtcΔt1Δt2+d2ΔtΔt1+dΔtΔt2=(Δt)22dΔtΔt1Δt2=(Δt)36}(10)

Solving Eq. (10) yields:

a=Δt3Δt12−4Δt3Δt1Δt2−6Δt2Δt12Δt2−4Δt3Δt22+12Δt2Δt1Δt22−24ΔtΔt12Δt22+24Δt13Δt2224Δt13Δt22b=Δt2Δt12−6Δt2Δt12Δt2+4Δt3Δt22−12Δt2Δt1Δt22+24ΔtΔt12Δt2224Δt13Δt22c=Δt3Δt1−2Δt3Δt2+6Δt2Δt1Δt212Δt12Δt22d=Δt26Δt1Δt2}(11)

Therefore, the proposed fractal scheme for Eq. (1) with second-order central space discretization is given as:

v¯i,jn+1=vi,jn+Δt1G(vi,jn,δxvi,jn,δyvi,jn,δy2vi,jn)(12)

v=i,jn+1=12(vi,jn+v¯i,jn+1)+Δt2G(v¯i,jn+1,δxv¯i,jn+1,δyv¯i,jn+1,δy2v¯i,jn+1)(13)

vi,jn+1=avi,jn+bv¯i,jn+1+cv=i,jn+1+ΔtdG(v=i,jn+1,δxv=i,jn+1,δyv=i,jn+1,δy2v=i,jn+1)(14)

where δxvi,jn=vi+1,jn−vi−1,jn2Δx,δyvi,jn=vi,j+1n−vi,j−1n2Δy,δy2vi,jn=vi,j+1n−2vi,jn+vi,j−1n(Δy)2.

Let G=β1∂v∂x+β2∂v∂y+β3∂2v∂y2 in Eq. (1), then the proposed scheme for this case can be written as:

v¯i,jn+1=vi,jn+Δt1(β1δxvi,jn+β2δyvi,jn+β3δy2vi,jn)(15)

v=i,jn+1=12(vi,jn+v¯i,jn+1)+Δt2(β1δxv¯i,jn+1+β2δyv¯i,jn+1+β3δy2v¯i,jn+1)(16)

vi,jn+1=avi,jn+bv¯i,jn+1+cv=i,jn+1+Δtd(β1δxv=i,jn+1+β2δyv=i,jn+1+β3δy2v=i,jn+1)(17)

3  Stability Analysis

The Fourier series analysis exists in the literature to find the stability conditions of partial differential equations. The scheme can be applied to non-linear PDEs, but it estimates the exact stability conditions in this case. For non-linear differential equations, equations are linearized first, and then this criterion is employed. This criterion is based on some transformations that convert the different equations into trigonometric equations with the involvement amplitude of waves. To apply this criterion, consider the following transformations:

v¯i,jn+1=P¯n+1eiIψ1ejIψ2,v¯i±1,jn+1=P¯n+1e(i±1)Iψ1ejIψ2v¯i,j±1n+1=P¯n+1eiIψ1e(j±1)Iψ2,vi,jn=PneiIψ1ejIψ2vi±1,jn=Pne(i±1)Iψ1ejIψ2,vi,j±1n=PneiIψ1e(j±1)Iψ2v=i,jn+1=P=n+1eiIψ1ejIψ2,v=i±1,jn+1=P=n+1e(i±1)Iψ1ejIψ2v=i,j±1n+1=P=n+1eiIψ1e(j±1)Iψ2}(18)

where I=−1. By substituting corresponding transformations from (18) into the first stage of the proposed scheme (15), it yields:

P¯n+1eiIψ1ejIψ2=PneiIψ1ejIψ2+Δt1{β1(e(i+1)Iψ1ejIψ2−e(i−1)Iψ1ejIψ22Δx)Pn+β2(eiIψ1e(j+1)Iψ2−eiIψ1e(j−1)Iψ22Δy)Pn+β3(eiIψ1ejIψ2−2eiIψ1ejIψ2+eiIψ1e(j−1)Iψ2(Δy)2)Pn}(19)

Simplifying Eq. (19) yields:

P¯n+1=Pn+Δt1{β1(eIψ1−e−Iψ12Δx)Pn+β2(eIψ2−e−Iψ22Δy)Pn+β3(eIψ2−2+e−Iψ2(Δy)2)Pn}(20)

Re-write Eq. (20) as:

P¯n+1=Pn+Δt1{β1Isinψ1Δx+β2Isinψ2Δy+β32(cosψ2−1)(Δy)2}Pn(21)

Let c¯1=β1Δt1Δx,c¯2=β2Δt1Δy,d¯1=β3Δt1(Δy)2. Then, Eq. (21) can be written as:

P¯n+1=Pn+{c¯1Isinψ1+c¯2Isinψ2+2d¯1(cosψ2−1)}Pn(22)

Similarly, substituting corresponding transformations from (18) into the second stage of the proposed scheme gives:

P=n+1=12(Pn+P¯n+1)+{c^1Isinψ1+c^2Isinψ2+2d^1(cosψ2−1)}P¯n+1(23)

By using Eq. (22) in Eq. (23), it gives:

P=n+1=12Pn+{12+(c^1Isinψ1+c^2Isinψ2+2d^1(cosψ2−1))}{1+c¯1Isinψ1+c¯2Isinψ2+ 2d¯1(cosψ2−1)}Pn(24)

Similarly, the following equation can be obtained by applying suitable transformations from (18) into the third stage of the proposed scheme (17):

Pn+1=aPn+bP¯n+1+cP=n+1+dΔt{c1Isinψ1+c2Isinψ2+2d1(cosψ2−1)}P=n+1(25)

where c1=β1ΔtΔx,c2=β2ΔtΔy,d1=β3Δt(Δy)2.

Upon substituting Eqs. (22) and (24) into Eq. (25) that gives:

Pn+1=Pn+b{1+c¯1Isinψ1+c¯2Isinψ2+2d¯1(cosψ2−1)}Pn+[c+dΔt(c1Isinψ1+c2Isinψ2+2d1(cosψ2−1))][12+{12+c^1Isinψ1+c^2Isinψ2+2d^1(cosψ2−1)}{1+c¯1Isinψ1+c¯2Isinψ2+2d¯1(cosψ2−1)}]Pn(26)

Eq. (26) can be written as:

Pn+1=(a¯+Ib¯)Pn(27)

where is a¯ and b¯ are real and imaginary parts of the expression in (26).

The amplification factor, in this case, is:

|Pn+1Pn|2=a¯2+b¯2≤1(28)

If the step size and βi′s,i=1,2,3 satisfy condition (28), the scheme will remain stable; otherwise, convergence is not guaranteed. Fig. 1 shows the amplification factor over the phase angle for both classical and fractal cases. The fractal scheme provides a smaller stability region for the chosen spatial step size than the scheme for the classical model. Also, for selecting small values of α, the stability region is becoming smaller and smaller.

Figure 1: Stability regions (a) classical (b) fractal models using β1=0=β2,β3=1,Δy=0.9, Δt(Δy)2=0.75

The next task is to provide convergence with the proposed scheme. For doing so, consider the system of fractal partial differential equations as follows:

∂u∂tα=A1∂u∂x+A2∂u∂y+A3∂2u∂y2+A4u(29)

where u is a vector and Ai′s,i=1,2,3 are matrices. Applying the proposed scheme on Eq. (29) as:

u¯i,jn+1=ui,jn+Δt1{A1δxui,jn+A2δyui,jn+A3δy2ui,jn+A4ui,jn}(30)

u=i,jn+1=12(ui,jn+u¯i,jn+1)+Δt2(A1δxu¯i,jn+1+A2δyu¯i,jn+1+A3δy2u¯i,jn+1+A4u¯i,jn+1)(31)

ui,jn+1=aui,jn+bu¯i,jn+1+cu=i,jn+1+Δtd(A1δxu=i,jn+1+A2δyu=i,jn+1+A3δy2u=i,jn+1+A4u=i,jn+1)(32)

Theorem 1: The proposed scheme (30)–(32) converges conditionally for Eq. (29).

Proof: Let the first stage of the exact scheme for (29) be:

U¯i,jn+1=Ui,jn+Δt1{A1δxUi,jn+A2δyUi,jn+A3δy2Ui,jn+A4Ui,jn}(33)

Subtracting Eq. (30) from Eq. (33) and considering:

U¯i,jn+1−u¯i,jn+1=e¯i,jn+1,Ui,jn−ui,jn=ei,jn,Ui±1,jn−ui±1,jn=ei±1,jn,Ui,j±1n−ui,j±1n=ei,j±1n

The resulting equation is:

e¯i,jn+1=ei,jn+Δt1{A1δxei,jn+A2δyei,jn+A3δy2ei,jn+A4ei,jn}(34)

Applying the norm on both sides of Eq. (34) and using triangle inequality yields:

e¯n+1≤en+c¯xen+c¯yen+4d¯yen+Δt1‖A4‖en(35)

where c¯x=‖A1‖Δt1Δx,c¯y=‖A2‖Δt1Δy,d¯y=‖A3‖Δt1(Δy)2.

The second stage of the scheme is expressed as:

u=i,jn+1=12(Ui,jn+U¯i,jn+1)+Δt2{A1δxU¯i,jn+1+A2δyU¯i,jn+1+A3δy2U¯i,jn+1+A4U¯i,jn+1}(36)

Similarly, the following inequality can be obtained by subtracting Eq. (31) from Eq. (36):

e=n+1≤12(en+e¯n+1)+c^xe¯n+1+c^ye¯n+1+4d^ye¯n+1+Δt2‖A4‖e¯n+1≤12en+(12+c^x+c^y+4d^y+Δt2‖A4‖)(1+c¯x+c¯y+4d¯y+Δt1‖A4‖)en(37)

Let the third stage of the exact scheme be:

Ui,jn+1=aUi,jn+bU¯i,jn+1+cU=i,jn+1+Δtd(A1δxU=i,jn+1+A2δyU=i,jn+1+A3δy2U=i,jn+1+A4U=i,jn+1)(38)

Similarly, subtracting the proposed scheme’s third stage from the exact scheme’s third stage and applying the norm yields inequality.

en+1≤aen+be¯n+1+ce=n+1+dcxe=n+1+dcye=n+1+d4dye=n+1+dΔt‖A4‖e=n+1(39)

≤aen+b(1+c¯x+c¯y+4d¯y+Δt1‖A4‖)en+(c+dcx+dcy+4dyd+dΔt‖A4‖)(12+(12+c^x+c^y+4d^y+Δt2‖A4‖)(1+c¯x+c¯y+4d¯y+Δt1‖A4‖))en+M(O((Δt)3,(Δx)2,(Δy)2))(40)

where cx=Δt‖A1‖Δx,cy=Δt‖A2‖Δy,dyx=Δt‖A3‖(Δy)2.

Re-write inequality (40) as:

en+1≤μen+M(O((Δt)3,(Δx)2,(Δy)2))(41)

where μ=a+b(1+c¯x+c¯y+4d¯y+Δt1‖A4‖)+(c+dcx+dcy+4dyd+dΔt‖A4‖)(12+(12+c^x + c^y+4d^y+Δt2‖A4‖)(1+c¯x+c¯y+4d¯y+Δt1‖A4‖)).

Let n=0 in inequality (41) it gives:

e1≤μe0+M(O((Δt)3,(Δx)2,(Δy)2))(42)

Since e0=0 because there is no error in initial conditions; therefore (42) is written as:

e1≤M(O((Δt)3,(Δx)2,(Δy)2))(43)

Suppose n=1 in inequality (41) then:

e2≤μe1+M(O((Δt)3,(Δx)2,(Δy)2))≤(1+μ)M(O((Δt)3,(Δx)2,(Δy)2))(44)

If this continued, then for finite n:

en≤(μn−1+…+μ+1)M(O((Δt)3,(Δx)2,(Δy)2))