A Finite Difference Method and Effective Modification of Gradient Descent Optimization Algorithm for MHD Fluid Flow Over a Linearly Stretching Surface

: Present contribution is concerned with the construction and application of a numerical method for the fluid flow problem over a linearly stretching surface with the modification of standard Gradient descent Algorithm to solve the resulted difference equation. The flow problem is constructed using continuity, and Navier Stoke equations and these PDEs are further converted into boundary value problem by applying suitable similarity transformations. A central finite difference method is proposed that gives third-order accuracy using three grid points. The stability conditions of the present proposed method using a Gauss-Seidel iterative procedure is found using Von-Neumann stability criteria and order of the finite difference method is proved by applying the Taylor series on the discretised equation. The comparison of the presently modified optimisation algorithm with the Gauss-Seidel iterative method and standard Newton’s method in optimisation is also made. It can be concluded that the presently modified optimisation Algorithm takes a few iterations to converge with a small value of the parameter contained in it compared with the standard descent algorithm that may take millions of iterations to converge. The present modification of the steepest descent method converges faster than Gauss-Seidel method and standard steepest descent method, and it may also overcome the deficiency of singular hessian arise in Newton’s method for some of the cases that may arise in optimisation problem(s).


Introduction
Magnetohydrodynamic (MHD) is also known as magneto fluid dynamics initiated by Nobel Laureate Hannes [Alfven (1942)] in 1970. Such fluids include electrolytes, plasmas, saltwater and liquid metals. The basic idea behind MHD is that currents can be induced in moving fluids by magnetic fields, which results in the polarisation of fluids and reciprocally alters the field itself. The application of MHD flows can be cited in the fields of Geophysics, Astrophysics, sensor Engineering and magnetic drug delivery. Numerous applications of MHD boundary layer flow over the surface in continuous motion has been reported in technical and industrial processes like in the field of metallurgy and chemical processes, moisture and thermal treatment of materials. For instance, in the polymer sheet extrusion from a die, in the extrusion of lamination and melt spinning process of polymers, fibres spinning, glass glowing and cooling of large metallic plates in a bath. Crane [Crane (1970)] was the first who studied the steady boundary layer laminar flows over a continuously stretching plate in which he solved the problem analytically and obtained an exact solution of Navier Stokes equations involved. Extension of the problem above was then made by Afzal et al. [Afzal and Varshney (1980)] to generalise the case with stretching velocity ~ where is the total distance from slit and is any constant parameter. Before Crane [Crane (1970)] who proposed the value of = 1 for velocity, Sakiadis [Sakiadis (1961)] limited the value of = 0 corresponding to the movement of the plate with constant velocity. Further development of such flows was made by Gupta et al. [Gupta and Gupta (1977)] on permeable stretching surfaces. They too obtained the exact solution of the flow phenomena and solution for thermal fluids were completely presented in Gamma function form. Ali [Ali (1995)] analysed the generalised form of such flows in which stretching velocity is taken in form of stretching sheets. Pavlov [Pavlov (1974)] highlighted the use of MHD application related to polymer processes and studied the boundary layer flow past a stretching plane sheet. Andersson [Andersson (1995)] extensively explained the similarity solution. He also studied that it is not only a similarity solution but an exact solution to the Navier Stoked equations completely. Extending the problem of Andersson et al. [Andersson (1995); Liu (2005)] studied the multiple effects for non-isothermal stretching sheet. He found the temperature distribution both in stated surface temperature and surface heat flux conditions, where n thermal condition and distance from the origin are linearly proportional. Due to a wide range of applications the unsteady boundary layer flow over a stretching sheet has been studied and analysed by Devi et al. [Devi, Takhar and Nath (1991); Elbashbeshy and Bazid (2004) ;Tsai, Huang and Huang, (2008)] have also studied the same problem above. In recent years, research has been measured for the investigation of two dimensional Burgers Huxley equation for understanding the various physical flows in fluid theory [Logan (1994); Granville and Burgen (2011)]. Thus, nowadays currently by implementing novel and a new approach to two dimensional Burgers Huxley equation yields in various finite difference, finite volume, open foam solver, Comsol Multiphysics, and finite elements approach [Taghizadeh, Akbari and Ghelichzadeh (2011);Ames (1965)]. Several powerful mathematical approaches such as Adomian decomposition method [Ismail, Raslan and Rabboh (2004) ;Hashim, Noorani and Batiha (2006)], spectral collocation method [Javidi (2006)], the tanh-coth method [Wazwaz (2008)], homotopy perturbation method [Hashemi, Daniali and Ganji (2007)], Exp-Function method [Zhou (2008)], variational iteration method [Batiha, Noorani and Hashim (2007)] and Differential Quadrature method [Sari and Güraslan (2009)] have been used in attempting to solve the Burgers-Huxley and the Huxley equations. Most of the solitary wave solutions of the generalized Burgers-Huxley equation have been studied by the learned researchers [Wang, Zhu and Lu (1990); El-Danaf (2007)]. Yi et al. [Yi and Chen (2012)] projected a Haar wavelet operational matrix method to solve a class of fractional partial differential equations. They develop the Haar wavelet operational matrix of fractional order integration and also fractional order differentiation. Using the operational matrix of fractional order differentiation, the fractional partial differential equations have been reduced to Sylvester equation. Wei et al. [Wei, Chen, Li et al. (2012)] extant a computational method for solving a class of space-time fractional convection-diffusion equations with variable coefficients which is based on the Haar wavelets operational matrix of fractional order differentiation. They also show error analysis in order to show the efficiency of the method. Previous study authors focused on spectral method and other numerical techniques for solving nonlinear and linear phenomena. It is worth mentioning that there is a vast amount of different approaches available in the literature to calculate the solutions of nonlinear systems of partial differential equations [Ray and Gupta (2013)]. Other numerical methods have also been employed, Haar wavelet method and variational iteration method were used [Macas (2017); Ray and Gupta (2013)]. Moreover, Burger fisher models and grey Scott models [Munafo (2011); Saqib, Hasnain and Mashat (2017)] investigated numerically ]. Recently, Shoaib et al. [Shoaib, Nawaz, Bibi et al. (2018)] study the flow of nanofluid because of heated disk rotation subjected to the convective boundaries with the chemical reaction of first order. The present attempt is made to propose a finite difference method to solve an MHD flow problem. The modification of an optimisation algorithm is presented and applied it to solve the difference equation obtained by applying the finite difference method on the transformed ordinary differential equation from the mathematical modelling of the fluid problem. The main advantage of the proposed finite difference method is that it is consisted on three grid points and provides accuracy of order three whereas the standard central difference formulas for first and second derivatives using three grid points provide the accuracy of order two. So the accuracy in order on fewer grid points can be enhanced by applying the presently proposed method for any order of numerical derivative. Besides, stability condition is found using Von Neumann stability criteria when the Gauss-Seidel method is applied for solving the difference equation. The leading error term arises in the Taylor series expansion of the obtained discretised equation is found, and the order of accuracy is also proved when the proposed finite difference method is applied to the ordinary differential equation. Immediate modification of the Gradient descent method is applied to solve the unconstrained optimisation problem constructed from the discretised equation. Present modification of the optimisation algorithm can be preferred to solve unconstrained optimization problem over some optimisation algorithms because it can provide quadratic rate of convergence due to the fact that it becomes Newton's method in optimisation when the value of the parameter contained in it is chosen to be zero, and for remaining non-zero values of the parameter, the algorithm may take more iterations to converge when it is compared with the Newton's method. However, this modified optimisation method gives more freedom to find the solution when the hessian becomes singular in Newton's method that may happen when Newton's method is applied to some problem in optimisation. So, the standard Newton's method has a deficiency to compute solution when the hessian becomes singular, but the present modification of the optimisation method may provide the solution in this/these situation(s). The comparison of the modified optimisation algorithm with Newton's method in optimisation and Gauss-Seidel iterative method is made by drawing tables. The maximum number of iterations, CPU time and norm of the error for two chosen values of the parameter contained in the modified optimisation algorithm are shown in the tables. These tables also show the maximum iterations, CPU time and norm of the error when the computations are performed using Newton's method in optimisation and Gauss-Seidel iterative method.

Problem formulation of fluid flow
Consider the flow over the flat plate, which is extended to infinity from its one end. The −axis is taken along the plate, and −axis is perpendicular on it. Using [Abdul (2013)], all derivatives in are neglected. It is further assumed that the flow is considered with the effect of the magnetic field 0 which is applied perpendicular to the − axis. Under the assumptions mentioned above, the governing equations of the flow can be expressed as: Subject to the boundary conditions = , = ( ), at = 0, > 0 is the velocity of the plate. In order to make the governing equations dimensionless, the following dimensionless variables are used [Abdul (2013)].
where 0 > 0 corresponds to the suction velocity of the plate and 0 < 0 corresponds to the injection velocity of the plate. By substitution the quantities (5)-(6) in Eqs.
(1)-(2) results in the following ordinary differential equation Since the term ′ in Eq. (7) is a dimensional term because ′ is not a constant, so in order to make the Eq. (7) dimensionless it is assumed that [26] ′ = , where is constant. By introducing into Eq. (7) results in the following differential equation ′′ + ( + 0 ) ′ − = 0 (8) The boundary conditions (3)-(4) of partial differential Eqs. (1)-(2) can be transformed to the following boundary condition of the ordinary differential equation

Numerical method
The boundary value problem (8) with boundary conditions (8a) is obtained from the considered flow problems (1)-(3). In order to solve this boundary value problems (8)-(8a) using the proposed numerical method, consider the standard fourth-order central difference formula for The formula (9) uses five grid points. If this formula is constructed on gird points −1 , , +1 then the central difference formula can be expressed by So, instead of considering the grid points −2 , −1 , , +1 , +2 on a fixed length of a subinterval ℎ, the half grid spacing is utilized in difference formula (10). But +1/2 and −1/2 are unknown quantities in the formula (10). So, in order to find these quantities, the interpolation of a parabola on grid points −1 , , +1 is taken into account. Let the equation of an interpolated shifted parabola is given by ≈ ( − ) 2 + ( − ) + (11) The equation of parabola (11) contains three unknowns. So, to find these unknowns, the Eq. (11) is evaluated on the points ( −1 , −1 ), ( , ) and ( +1 , +1 ). By substituting these mentioned three points in Eq. (11) gives the three equations in three unknowns and solving these equations gives the values of the unknown , and which are expressed as By evaluating the Eq. (11) using the values of , and mentioned above and introducing −1/2 = + +1 2 and +1/2 = + +1 2 into the resulting equation gives Since the formulas (12a) and (12b) are third-order accurate, so the resulting central differencing formula (10) using formulas (12a)-(12b) is also third-order accurate.

Stability of the proposed method
The stability of the proposed method using a Gauss-Seidel iterative procedure is checked by employing Von-Neuman stability criteria on Eq. (16) using Eqs. 17(a)-17(b). In order to begin the stability procedure of Von-Neuman stability criteria, the following substitutions are considered Under these substitutions, Eq. (16) with the substitution of Eqs. 17(a)-17(b) yields the following equation and considering this stability procedure for the special case when = 0 Rearranging and dividing the resulting equation by yields According to Von-Neuman stability criteria, the resulting Eq. (18) must satisfy By using the stability condition (19) for Eq. (18), the resulting equation is given as: So, when = 0 and the parameters ℎ, 0 and satisfy the inequality (21), the presently proposed methods (16)-(17b), using Gauss-Seidel iterative procedure remains stable.

Error analysis
In order to find the leading error term and order of the proposed method, the following procedure is considered. Rewrite the Eq. (16) in the following form Since +1/2 + −1/2 ≈ 1 8 By substituting the results from Eqs. (25)-(26) into Eq. (22) yields Since ′′ + 0 ′ − = 0. So Eq. (30) can be expressed as: The leading error term can be expressed by

Modified optimization algorithm
The standard Gradient descent method can be used to solve the difference (Eqs. (16)-17(b)) obtained by applying the proposed third-order central difference method on the boundary value problem obtained from the considered flow problem. However, due to its low rate of convergence, the solution can be obtained with a large number of iterations which may converge using a small fixed step size. One of the advantages in using the Gradient descent method is that it is simple and easy in implementation for the problem when compared with the Conjugate Gradient Algorithm (CGA), but the CGA converges faster than the Gradient Descent Algorithm (GDA). The deficiency in applying the standard GDA is its lower rate of convergence, and this results in the consumption of time to get the solution of the problem with the fixed step size. In order to increase the rate of convergence of GDA, the standard GDA can be improved, and the resulted modified optimisation algorithm is shown further in the present work provide convergence near to Newton's method. Newton's method in the optimisation uses the hessian, and the problem of finding the solution using Newton's method in optimisation is that the hessian may become singular that may happen in some of the case(s). However, the presently modified approach of the optimisation algorithm uses the term consisted of the sum of a constant multiple of an identity matrix and the hessian. So the singularity of hessian may have been controlled by the parameter contained in the modified descent algorithm. The modified descent gradient algorithm in optimisation is constructed in the following Theorem. Theorem. The iterative formula for the presented modified gradient descent algorithm (MGDA), is given by +1 = − ( + )∇ ( ) where = ∇ 2 ( ) is the Hessian, ∇ ( ) denotes the Gradient and is the step size. Proof. Consider the Euclidean manifold and the squared length of small incremental vector that connects two points and on this manifold, and this can be expressed as In order to find the gradient descent direction, consider the direction that minimizes ( + ) under the constraint ‖ ‖ 2 = 2 , (33) Which is defined by the vector . In (33) is very small constant. Let = (34) The task is to find a vector that minimizes under the constraint ‖ ‖ 2 = ∑ 2 =1 = 1.
(36) The standard descent method uses the first-order Taylor approximation of ( + ), but the present modified form of gradient descent method consists of the extended form of Taylor series for ( + ) Given in (35). By using the optimal approach, the Lagrange method is utilised to find an optimal direction . In order to find optimal direction , the following equation is constructed By applying partial derivative on Eq. (37) gives where is a Lagrange multiplier. By taking the vectored derivative in (38) and solving the resulting equation for the direction yields = −[ 2 ( ) + ] −1 ( ) (39) By using (39), the iterative formula for the modified gradient descent method can be expressed as: In order to implement the presently constructed modified gradient descent algorithm on the difference (Eqs. (12a)-(12b), (15)), the following metric on the Euclidean manifold is defined.
By substituting +1/2 and −1/2 from (12a)-(12b) into Eq. (41) yields the following metric on Euclidean manifold where is the objective function. The gradient or derivative of the metric defined in (42) is given by denotes the Jacobian, and it is defined as: Since the presently modified Gradient descent method used the hessian in its every iteration, so the hessian is given by

Results and discussions
This work presents the third-order central finite difference method for the flow over the flat plate with the effect of the magnetic field. The order of accuracy can be verified from the error analysis mentioned above that yields a leading error term containing the product of ℎ 4 . The consistency of the proposed method can be verified by applying the limit as ℎ → 0 to Taylor series expansion form of the difference equation; if the original equation is obtained, the method will be stable otherwise it will be unstable. So, it can be verified that the proposed finite difference method is consistent. Also, since the method is more than firstorder accurate, so it is consistent. Therefore, the stability condition (21) can be considered as the conditions for convergence of the method. Since the method is conditionally stable and consistent, so it is conditionally convergent for the considered constructed flow problem. The attained stability condition was founded when the Gauss-Seidel iterative procedure applied to the discretised form of the problem (8) using the third-order central proposed method. Thus, the above-founded stability condition (21) is valid for the case of Gauss-Seidel iterative process. However, the consistency criteria are fixed for any iterative procedure because it depends upon the order of the method and the method is third-order accurate regardless of any iterative process because the order of the method cannot be affected by an iterative method. Also, it is to be noted that the number of iterations of the present method using a Gauss-Seidel iterative method can be reduced by using the Newton iterative method for difference (Eqs. (16)-(17b)). Newton's method takes only a few iterations to produce the required results, but the standard Gauss-Seidel iterative procedure took hundreds or thousands of iterations to converge with the given stopping criteria in some of the cases. The stopping criteria for the iterative method are given by The solutions of boundary layer flows must have asymptotic behaviour at one end. Moreover, since the problem is constructed on the semi-infinite domain, which is restricted to some finite part of the domain provided that the obtained solution must asymptotically approach to zero. It has been observed in most of the present cases of the simulations that the solution converges over the large transformed domain ( ) and it can satisfy the asymptotic behaviour characteristic over the large values of the independent variable. Figs. 1-3 are drawn for velocity profile with the variations of parameters given in Eq. (8) also, the norm of error over the log of the number of iterations. Since it has been given in Abdul [Abdul (2013)] that the steady flow corresponds to the case when = 0 and the case ≠ 0 represents the unsteady flow. Fig. 1 shows the variations of velocity profile with the increasing values of suction/injection parameter 0 . The increasing values of the parameter 0 affect the boundary layer's thickness. The thickness of the boundary layer increases by enhancing the values of the parameter 0 for both suction and injection case but maintains an exponentially decaying behaviour. Fig. 2 shows the velocity profile with the variation of the magnetic field parameter. Since enhancing the applied magnetic field on electrically conducting fluid generates the Lorentz force that results in a reduction in the velocity of the fluid. The thickness of the boundary layer also decreases with the increase of the magnetic parameter, and this can be seen in Fig. 2. Fig. 3 depicts that the thickness of the boundary layer reduces with the increase in parameter . For steady flow, the thickness of the boundary layer is larger than the thickness of unsteady flow and also velocity profile decreases when the parameter is enhanced. Fig. 4 shows the skin friction coefficient with the variation of suction/injection velocity parameter. Since, the slope of the tangent at point zero on velocity profile decreases with the increase in the suction/injection velocity parameter due to decrease in the velocity profile, so negative of the tangent's slope at zero increases. The increasing behaviour of the skin friction coefficient can be seen in Fig. 4. Similarly, using the slope of the tangent concept from calculus on the velocity profile at zero also have similar behaviour that has been discussed earlier when the magnetic parameter and parameter varies and the increase in the skin friction coefficient with the variation of magnetic parameter and the parameter can be seen in Figs. 5-6. In order to find the derivative numerically at initial point zero, forward sixth-order difference formula is implemented at point zero. Fig. 7 shows the norm of the difference of the numerical solution obtained from the presently proposed third-order central difference method and optimisation algorithm modified in the present work over the log of the number of iteration with the variation of the number of nodes. When the number of nodes increases, the maximum number of iterations taken by the finite difference using the modified optimisation algorithm also increases with the stopping criteria. However, it should be noted that the maximum number of iterations also depends upon the parameter contained in the optimisation algorithm modified in the present contribution. If this parameter decreases and approaches to zero, the maximum number of iterations will have non-decreasing behaviour. Present modification of optimisation algorithm can give approximately the same maximum number of iterations which are obtained from Newton's method. Fig. 8 shows the maximum number of iterations taken by the finite difference method using modified steepest descent method over the step size. It can be seen that the maximum number of iterations increases with the increase in the step size. Tab. 1 shows the norm of the error between the exact solution and the numerical solution obtained from the finite difference optimization algorithm using three methods. It can be seen from Tab. 1 that this norm of the error is the same as the one obtained by three methods. Tab. 2 shows the comparison of three algorithms in terms of the maximum number of iterations and CPU time consumption. It can be seen from Tab. 2 that the present modification of the steepest descent method converges faster than the Gauss-Seidel iterative method but little expensive than Newton's method.
The two advantages of the present modification of the optimization algorithm are that it can converge with the rate of convergence of order two and also it can provide more freedom to find the solution when the Hessian becomes singular. These are the main reasons when this optimization method can be considered as the replacement of Newton's method and some other algorithms due to its faster convergence and the problem of the singularity of hessian.

CPU performance
Since the maximum number of iterations are just two for the case of the modified optimization algorithm, and the time spent is also shown in Tab. 2. Due to these facts, the CPU performance is not much affected by similar kind of simulations. However, the CPU time may vary when the simulations are performed using a higher number of nodes with some small value of the parameter near to zero but not very small. Since the laptop HP CORE i5, 2.4 GHz processor, 2 Cores, 4 GB RAM is used during the simulation which has some heated problem during high simulations (those simulations which may take large time) to get the results in some or most of the time taking simulations. However, those results are not displayed in the present contribution. The elapsed time is shown in Tab. 2 is just fraction of a second taken by the modified optimization algorithm and Newton's method to get the result under the given circumstances and this elapsed time is found from the code constructed in Matlab. In addition to the time shown in Tab. 2, some additional information can be obtained using Matlab when the code is run in Matlab.
Since the little higher value of the parameter contained in the modified optimization algorithm leads to the higher number of iterations, so the time and memory can be checked, and it is displayed in Tab. 3. Step size Max. no. of iters.

Conclusion
A finite difference method with stability and error analysis has been presented for solving the electrically conducted fluid flow problem over the linearly stretching sheet with the effect of the magnetic field. An optimisation algorithm that is extended form of standard descent method is proposed and applied this extended algorithm to solve the difference equation obtained from the fluid problem when it is discretised by the presently proposed central difference method. It has been shown that the modified optimisation algorithm is faster than the Gauss-Seidel method to solve the difference equation, and both methods give the same norm of the error. The present modification of the standard descent method can give the quadratic rate of convergence because it becomes Newton's method when the value of the parameter contained in it is chosen to be zero and also it may overcome the deficiency of Newton's method when the Hessian becomes singular in some of the case(s). Under the advantages mentioned earlier, this present modified extension of standard descent method can be applied further to solve optimisation problems, and this method can be the excellent replacement of Newton's method and others optimisation algorithms. The present finite difference method can be similarly extended to higherorder when the interpolation of any higher-degree polynomial can be considered to find the unknown points evaluated on half grid spacing.