1 Introduction

The Jeffery-Hamel problem is related with the incompressible viscous fluid flows between non-parallel walls. It has been broadly utilized in the applications of fluid mechanics,e.g.,civil, environmental,mechanical,and bio-mechanical engineering. The origin and mathematical description for such a type of fluids has been introduced by Jeffery^[1] and Hamel^[2]. The Jeffery-Hamel flow presents an exact similarity solution of the Navier-Stokes equations in the special case of two-dimensional (2D) flows through the channels with inclined plane walls,meeting at the vertexes and with source or sink at the vertexes. This issue has been widely studied by many authors^{[3, 4, 5, 6, 7]}. The classical Jeffery-Hamel problem was extended in Ref.^[8] to study the effects of the external magnetic field on conducting fluids. The magnetic field works as a control parameter,along with the Reynolds number of the flow and the angle of the walls.

To obtain the exact solution in the magnetohydrodynamics (MHD) Jeffery-Hamel problem,one has to deal with a dynamical system with strong nonlinearity. Many different analytical and numerical methods have been proposed to solve such complex problems,e.g.,the Adomian decomposition method^{[9, 10]},the homotopy perturbation method (HPM)^[11],the differential transform method (DTM)^[12],the homotopy analysis methods (HAM)^[13],and the variational iterational method (VIM)^{[14, 15]}. However,the approximate solutions for the Jeffery-Hamel flow problems have seldom been reported^{[16, 17, 18, 19]}.

The strength of the artificial intelligence algorithms based on the modern deterministic procedures of the supervised learning and unsupervised artificial neural networks (ANNs) optimized with the stochastic techniques and the swarm intelligence algorithms has been exploited in a variety of applications^{[20, 21, 22]}. Recently,these techniques were used to solve nonlinear oscillators^[23],Lane Emden flower equations^[24], one-dimensional (1D) nonlinear Bratu's problems^{[25, 26]},boundary value problems (BVPs) of 2D nonlinear Bratu equations^{[27, 28, 29]},nonlinear equations^[30],fluid dynamic problems^{[31, 32]}, mathematical models of thin film flows of the third-grade fluids^[33],nonlinear BVPs of Troesch's equations^{[34, 35]},BVPs of pantograph functional differential equations^[36],initial value problems of nonlinear Painlev\'{e} equations^{[37, 38]},and linear and nonlinear fractional differential equations^{[39, 40]}. Besides, the solution of the fractional order system based on the Riccati and Bagley-Trovik fractional differential equations is another illustrative application of such solvers^{[41, 42, 43]}. All these are the motivation factor for the authors to investigate such techniques and to provide an alternate,accurate,and effective platform for solving the variants of the MHD Jeffery-Hamel problems.

In the present study,a new stochastic technique is developed based on the neural networks optimized with the evolutionary calculation and the sequential quadratic programming techniques for the solution of nonlinear Jeffery-Hamel flow equations. The original 2D problem is transformed into an equivalent higher order BVP of ordinary differential equations (ODEs). The mathematical model of the transformed equation is constructed with the help of feed-forward ANNs in an unsupervised manner. The weights of the networks are obtained by the evolutionary computation based on the genetic algorithms (GAs) and the SQP method. The designed scheme is evaluated on the variants of the Jeffery-Hamel problem. The proposed solutions are compared with the numerical results of MATHEMATICA obtained by the fully explicit Runge-Kutta method,the DTM^[44], HPM^{[44, 45, 46]},the HAM^{[44, 45, 46, 47]},the VIM^[46],and the optimal homotopy asymptotic method (OHAM)^[48]. Moreover,a sufficient large number of independent runs of the proposed schemes are executed to validate the accuracy, convergence,and effectiveness of the proposed algorithm.

2 MHD Jeffery-Hamel problem

In this section,a brief description for the transformation of a 2D MHD Jeffery-Hamel flow problem based on partial differential equations into third-order BVP ODEs is presented. Many authors have reported the solution for the transformed problem with the help of different numerical and analytical techniques,which are always simpler to handle than that of the original partial differential equations^{[6, 44, 45, 46, 47, 48]}.

2.1 Governing equations

The equation governing the flow of an incompressible MHD fluid and neglecting the gravity and thermal effects is

where ${V}$ is the velocity vector defined by

$V=(u(x,y,z),v(x,y,z),w(x,y,z)),$

$\rho $ is the constant density, ${\rm D}/{\rm D}t$ denotes the material time derivative,\textit{p} is the pressure,$T$ is the extra stress tensor,$J$ is the current density,and $B$ is the total magnetic field.

For the incompressible Newtonian fluid,

where $\mu $ is the constant viscosity,and ${A}_{1 }$ is the first Rivlia-Ericken tensor given by

2.2 Formulation of problem

Consider a 2D steady and fully developed laminar flow of an incompressible conducting Newtonian fluid from the source/sink between the two plates meeting at the angle 2$\alpha$ (see Fig. 1). The flow can be either converging or diverging. As shown in Fig. 1, the cylindrical polar coordinate $(r,\theta,z)$ is considered,and the two plates intersect in the $z$-axis.

The flow is assumed to be purely radial with no slip condition at $\theta = \pm \alpha$ and $v = w = 0$. Then,the seeking velocity profile is given by

Assuming that the magnetic field $B$ is perpendicular to the velocity field $V$,and the induced magnetic field is negligible as compared with the imposed magnetic field. Then,the momentum equation given in (2) can be reduced to

where $\nu$ is the coefficient of the kinematic viscosity,i.e., $$\nu =\frac{\mu}{\rho},$$ $u(r,\theta$) is the velocity along the radial direction,$\sigma$ is the conductivity of the fluid,and $B_{0}$ is the electromagnetic induction.

The continuity equation given in (1) can be reduced to

which requires $ru(r,\theta)$ to be independent of $r$. Substituting into (7) gives

This implies where $c$ is a function of $r$ only. Substituting (9)--(10) in (6) and differentiating the obtained results with respect to $\theta$ yield

Introducing the dimensionless parameters in (11),i.e.,

yields the transformed problem as follows: where $Re$ and $Ha$ are the Reynolds number and the Hartmann number,respectively,which can be defined by

In the above equations,$u_{\rm {max} }$ is the velocity at the center of the channel ($r = 0$),and

3 Neural network mathematical modeling

Neural networks are well-known for their competency,used as a universal function approximation. In the ANN methodology,the solution $y(\eta)$ and its $n$th-order derivative ${\rm d}^{n}y(\eta)/({\rm d}t^{n})$ of the differential equation are approximated with the following continuous mapping^{[49, 50]}:

where $\delta_{i}$,$w_{i}$,and $\beta_{i}$ are the $i$th elements of the bounded real-valued adaptive weights,and $m$ is the number of the neurons in the ANN architecture. $f$ is the activation function,and it is normally taken as the log-sigmoid function expressed as

for the hidden layers and the linear function for the output layer. With the log-sigmoid function, we can obtain the ANNs for the solution $y(\eta)$ with the first-order derivative (${\rm d}y(\eta)/({\rm d}\eta)$) and the third-order derivative (${\rm d}^{3}y(\eta)/({\rm d}\eta^{3})$) as follows:

The arbitrary combinations of the networks given in (15)-(17) are used to develop an approximate model of the problem (12) and its architecture. The results are shown in Fig. 2.

The fitness function $\epsilon$ is formulated by defining an unsupervised error based on the sum of two mean square errors as follows:

The error term $\epsilon_{1}$ for (12) is

where $\eta\in(0,N)$ is divided into $K$ intervals,i.e., $\eta\in(\eta_0,\eta_1,\cdots,\eta_K)$ in which $\eta_0=0$ and $\eta_K=N$,and

Similarly,the error term $\epsilon_{2}$ for the initial and boundary conditions for (12) can be given by

It is quite clear that the error terms $\epsilon_{1}$ and $\epsilon_{2}$ approach zero with the proper parameters $\delta_{i}$,$w_{i}$,and $\beta_{i}$,and thus the unsupervised error $\epsilon$ approaches zero. Consequently,the exact solution $y(\eta)$ is closely approximated with the proposed solution $\hat{y}(\eta)$.

4 Learning methodologies

In this section,a brief description of the methodologies used for the weight training of neural networks is presented,including the GA,SQP,and GA-SQP algorithms. These algorithms are described with the help of the generic flow diagrams,procedure steps,and parameter settings.

The SQP method belongs to a class of nonlinear programming techniques,used excellently for constraint optimization problems. Its supremacy is well recognized in terms of the efficiency, accuracy,and percentage of the successful solutions over a large number of test problems. An extensive introduction and an overview of the SQP methods have been given by Nocedal and Wright^[51]. Necessary detail and mathematical background about the SQP methods can also be found in Refs.^[52] and ^[53]. The SQP methods have been used extensively in diverse fields of engineering and applied sciences. However,few significant applications have been addressed with the SQP algorithm^{[54, 55]}.

The evolutionary computing techniques mainly based on the GAs are inspired by the natural evolution process. The GAs were introduced by Holland in 1975 in his pioneer work to mimic a simple picture of natural selection^[56]. The performance of the GAs depends upon the diversity of the initial population,the suitable selection of the fittest chromosomes to the next generation,the survival of good genes in the recombination operation,and the seeding of new genetic materials in mutation. The GA based on the evolutionary algorithms is broadly used due to its capacity of controlling the robustness and avoiding the local convergence,applicability in divergent environment,and higher efficiency and effectiveness compared with other meta-heuristic solvers^{[57, 58, 59, 60]}.

Besides the GA and SQP methods,the hybrid approach GA-SQP is also used for the weight training of the neural networks to optimize the fitness function $\epsilon$ given in (18). The MATLAB of version 2011a is incorporated for running the routines for the GA and SQP methods. The generic flow diagram of the proposed hybrid scheme GA-SQP is given in Fig.3(a),and the usual working reproduction operators of the GA are given in Fig.3(b).

In order to find the desired weights of the neural networks,a stepwise procedure for the GA-SQP approach is given as follows:

Step 1 Initialization

Generate the bounded real values to create the initial population that consists of a number of chromosomes or individuals. Each chromosome or individual has a number of genes,which is equal to the number of the weights in the neural network model of the equation. A population of $N$ chromosomes is mathematically given by

and the $i$th chromosome ${C_{{r},i}}$ is presented mathematically by

where $m$ is the number of the neurons of the neural network model. The necessary parameter settings used for the execution of the GAs are given in Table1.

Step 2 Fitness evaluation

Calculate the fitness values for each individual of the population by use of the expression given in (18) for the fitness function.

Step 3 Termination criteria

Terminate the iterative process of the algorithm on the basis of the criteria that the value of the fitness function $\epsilon$ is reduced to the predefined value,and the predefined number of the generations is exceeded. If the termination criterion is fulfilled, then go to Step 6.

Step 4 Ranking

Rank the individuals of the populations on the basis of $\epsilon$.

Step 5 Reproduction

Reproduce the population for the next generation by applying the genetics search operators such as the crossover, mutation,selection,and elitism. A general working of reproduction operators is given in Fig.3(b).

Go to Step 2.

Step 6 Refinement

Use the SQP local search algorithm for the speedy optimization by taking the best individual of the GAs as the starting point or the initial weights of the algorithm. The optimization parameters for the SQP are set as those listed in Table1.

Step 7 Storage

Store the best individual,i.e.,the optimal weights obtained for this run of the GA-SQP algorithm. Repeat the procedure from Step 2 to Step 6 for a sufficiently large number of times for the reliable statistical analysis.

5 Simulations and results

In this section,the results of the proposed scheme are presented for the two nonlinear Jeffery-Hamel problems. Five variants of each problem are investigated by taking different values of the Reynolds number,the Hartmann number,and the angle $\alpha$. The comparative studies of the proposed results are tabulated with the reference numerical and reported solutions for each case.

5.1 Problem I: variants of Jeffery-Hamel flow problem based on ${{Re}}$ and $\alpha$

The governing equation for the Jeffery-Hamel problem for this case is obtained from (12) by setting $H = 0$,i.e.,

The designed scheme is evaluated on five test cases of the problem, i.e.,Case I^[44]:

Case II^[44]:

Case III^[45]:

Case IV^[48]:

and Case V: a variety of values of $Re$ and $\alpha$. The exact solution of the problem is not known for each case. Therefore,the numerical solutions have determined for each case of the problem by use of the inbuilt routine of MATHEMATICA based on the fully explicit Runge-Kutta method,and these solutions are used as the standard or reference for comparison.

We have used the neural networks with 10 numbers of neurons to model the problem. Therefore,a total of 30 unknown parameters,i.e., $(\delta_{i},w_{i},\beta_{i}$),are required for the optimization of the networks. The fitness function $\epsilon$ is formulated for (21) by taking the inputs between $\eta \in ^{[0, 1]}$ with a step of 0.1,i.e.,$K = 10$,as follows:

The proposed schemes based on the GA,SQP,and GA-SQP are applied to optimize (22) with the parameter settings given in Table1. One set of weights for the respective GA,SQP,and GA-SQP algorithms with the fitness values $\epsilon$ of

for Case I,

for Case II,

for Case III,and

for Case IV are given in Tables A1-A4 in Appendix A. These weights can be used in (15) for calculating the solution of the problem for any inputs between 0 and 1. The solutions determined along with the values of the absolute error ($A_{\rm e}$) from the reference results for the inputs $\eta\in [0,1]$ with a step of 0.1 are provided in Tables 2-6 for four variants of the Jeffery-Hamel flow by the proposed schemes.

The reported solutions of the three analytical methods,i.e.,the DTM^[44],the HPM^[44],and the HAM^[44],are also given in Tables 2--5 for Case I and Case II. The reported results of the HPM^[45] for Case III are provided in Table6. The quoted solutions with their absolute errors for the VIM^[46], OHAM^[48],and HPM^[46] are provided in Table7. The mean absolute errors ($G_{\rm {MAE}}$) for the DTM,HPM,HAM,GA,SQP, and GA-SQP algorithms are

respectively,for Case I and

respectively,for Case II. The mean absolute errors for the HPM,GA,SQP,and GA-SQP algorithms are

respectively,for Case III. The mean absolute errors are

for the reported solutions of the OHAM,HPM,and VIM,respectively,while are

for the proposed GA,SQP,and GA-SQP techniques,respectively,for Case IV. In general,the results of the proposed scheme based on the SQP and GA-SQP techniques show good agreement with the standard numerical solutions,and are relatively better for the GA-SQP approach.

For Case V,the variants of the problem (22) are solved by the proposed hybrid method GA-SQP by taking one parameter of $Re$ or $\alpha$ fixed and varying the other parameters. The proposed velocity profiles are shown in Fig.4 for the inputs between [0,1] with a step size 0.05 along with the reference numerical solution for each case. Moreover,the achieved values of the mean absolute error and fitness are plotted in Fig.5 for each variant. It is seen that the proposed results are consistently overlapping with the reference solution with good accuracy.

5.2 Problem II: variants of Jeffery-Hamel flow problem based on ${{Ha}}$, ${Re}$,and $\alpha$

The governing equation for this case study is

We will investigate for the solution of the problem for four variants with $Re = 50$,i.e.,Case I: $Ha = 250$ and $\alpha= 7.5^{\circ}$^[45],Case II: $Ha = 500$ and $\alpha= 7.5^{\circ}$^[45],Case III: $Ha = 1\,000$ and $\alpha = 5^{\circ}$^[47],and Case IV: $Ha = 1\,000$ and $\alpha = -5^{\circ}$^[47]. For Case V,various values of $Re$,$Ha$,and $\alpha$ are considered. The reference numerical solutions are calculated for each variant of the problem by use of the Runge-Kutta method.

We have applied the designed methodology on a similar procedure as in the previous problem. However,the fitness function $\epsilon$ for the inputs of the training set between $\eta\in [0,1]$ with a step of 0.1 is formulated for (23) and $Re= 50$ as follows:

Now,the requirement is to find the weights for the optimization of the fitness function $\epsilon$ given in (23) by the GA,SQP,and GA-SQP algorithms with the values of the parameters given in Table 1. A set of the trained weights for the respective GA,SQP,and GA-SQP methods with the fitness values $\epsilon$ of 1.7427$\times10^{-2}$,1.6305$\times10^{-9}$,and 2.1525$\times10^{-9}$ for Case I,1.1500$\times10^{-2}$, 3.2118$\times10^{-9}$,and 1.7749$\times10^{-8}$ for Case II, 1.5915$\times10^{-2}$,2.0590$\times10^{-8}$,and 1.9414$\times10^{-8}$ for Case III,and 6.0675$\times10^{-3}$, 1.8927$\times10^{-9}$,and 1.7381$\times10^{-8}$ for Case IV are given in Tables A5-A8 in Appendix A. These weights are used to calculate the solutions of the four cases of the problems,and the results are tabulated,respectively,in Tables 8--11 along with the values of the absolute errors for the inputs $\eta\in ^{[0, 1]}$ with a step of 0.1.

The reported solutions along with the values of $A_{\rm e}$ for the analytical technique HPM^[45] are also given in Tables 8 and 9 for Cases I and II,respectively. The reported results of the HAM^[47] for Cases III and IV are provided in Tables 10 and 11, respectively,for the inputs $\eta\in^{[0, 1]}$ with a step of 0.1. The values of $G_{\rm {MAE}}$ for the HPM,GA,SQP,and GA-SQP methods are

respectively,for Case I,

respectively,for Case II,

respectively,for Case III,and

respectively,for Case IV. In general,the results of the SQP and GA-SQP algorithms are found to be closely matched with the standard numerical solution.

For Case V,the variants of the problem (23) are solved by the proposed hybrid method GA-SQP by fixing any two parameters of $H$, $Re$,and $\alpha$ and varying the third one for both convergent and divergent fluid flows. The proposed velocity profiles are shown in Fig.6 for the inputs between [0,1] with a step size 0.05 along with the reference numerical solution for each variant of the problem (23). Moreover,the values of the mean absolute error and fitness are plotted in Fig.7 to observe the level of accuracy. It is seen that the proposed results are consistently matching with the reference solution with seven to eight decimal places of accuracy.

6 Statistical analysis for solvers

The reliable and effective inferences for the results of the stochastic solvers can only be made on the basis of a large number of independent runs rather than one single successful execution of the algorithm. In this section,the results are presented for the proposed schemes based on a sufficiently large number of independent runs of the algorithms to examine the accuracy,convergence,and computational complexity.

The accuracy of the solvers is analyzed based on the statistical parameters,i.e.,the mean absolute error ($G_{\rm {MAE}}$),the standard deviation of the absolute error ($G_{\rm {STD}}$),and the minimum absolute error ($G_{\rm {MIN}}$) from the standard numerical results. The results are calculated for the SQP and GA-SQP methods based on 100 independent runs for the inputs between 0.1 and 0.9 with a step of 0.2 for each case of both problems (see Table12). It is observed that there is no significant difference in $G_{\rm {MIN}}$ for the SQP and GA-SQP algorithms for different variants of the MHD Jeffery-Hamel flow. However,the values of $G_{\rm {MAE}}$ and $G_{\rm {STD}}$ for the hybrid GA-SQP technique are invariably much better than those for the SQP algorithm.

The values of the fitness $\epsilon$ and $G_{\rm {MAE}}$ for 100 independent runs of the GA,SQP,GA-SQP algorithms are calculated. The results of $\epsilon$ and $G_{\rm {MAE}}$ are plotted against the number of the independent runs of the algorithms in Fig.8(a) and Fig.8(b),respectively,for Case I to elaborate the small difference in the value results on the semi-log scale. It is observed that with the increase and decrease in the fitness,$G_{\rm {MAE}}$ varies accordingly. The values of $G_{\rm {MAE}}$ are plotted against the number of the independent runs for different cases of Problems I and II in Fig.8(b) to Fig.8(i), respectively. Few effects can be made from the results. First,the hybrid GA-SQP algorithm provides consistently lower values of $G_{\rm {MAE}}$ than those of the GA and SQP techniques. Secondly, for few iterations of the SQP algorithm,the divergence is observed almost in all cases of the problem. It is quite understandable because the SQP algorithm is a local optimization algorithm,for which the guarantee cannot be given for the consistent convergent results. Thirdly,the results of $\epsilon$ and $G_{\rm {MAE}}$ of the GAs are relatively much more inferior than those of the SQP and GA-SQP techniques.

The reliability of the proposed solvers is examined on the basis of the percentage of the convergent runs,i.e.,given by the pre-defined values of $\epsilon$ and $G_{\rm {MAE}}$. The results of the convergence analysis based on 100 independent runs of the GA,SQP and GA-SQP algorithms are provided in Table12. On the basis of the fitness $\epsilon\leqslant 10^{-6}$,the average convergence rates for the SQP and GA-SQP algorithms are 96.75% and 98.75% for Problem I and 91.25% and 94.50% for Problem II. Whereas based on $G_{\rm {MAE}}$ less than $10^{-6}$,the average convergence rates for the SQP and GA-SQP algorithms are 78.75% and 85.25% for Problem I and 65.50% and 66.00% for Problem II.

We have analyzed the proposed solvers further based on the mean absolute error (G$_{\rm MAE}$) and the mean fitness ($M_{\rm F}$) defined by

where $P$ and $R$ are integers,representing the total number of the grid points and the total number of the independent runs of the solver,respectively. $y_{i}$ denotes the numerical solution for the $i$th input,$\hat{y}_{i}^{r}$ is the proposed approximate solution for the $i$th input and $r$th independent run,and $\epsilon_{{ r}}$ represents the fitness value for the $r$th run of the algorithm. In our experimentation,we have taken $\eta\in[0,1]$ with a step of 0.1,i.e.,$$P = 11,$$ and 100 independent runs,i.e.,$$R = 100.$$ The values of $G_{\rm MAE}$ and $M_{\rm F}$ along with $G_{\rm {STD}}$ are determined for the GA,SQP,and GA-SQP algorithms. The results are presented in Table13 for each case of the problems. The values of $G_{\rm MAE}$ and $M_{\rm F}$ are found to be the best for the results due to the GA-SQP technique than those of the GA and SQP schemes for all cases of the BVPs.

The computational complexity of the schemes is analyzed based on the computational time taken for finding the optimal weights of the neural networks by each algorithm. The analysis is performed based on 100 independent runs of each technique. The results are also provided in terms of the mean execution time ($t_{\rm {MET}}$) along with the standard deviation of the absolute errors ($G_{\rm {STD}}$) listed in Table14. It is found that the values of $t_{\rm {MET}}$ for the GA-SQP algorithm are larger than those of the GA and SQP algorithms. However,this aspect can be compromised based on the supremacy of the GA-SQP algorithm in terms of the convergence and accuracy over other two techniques. The results of the computational time analysis are based on a Dell Precision 390 Workstation with the Intel (R) Core (TM) 2 CPU 60002.40GHz processor. The program is running with the MATLAB of version 2011a.

7 Conclusions

On the basis of the results presented in the last sections,the following conclusions can be made:

(i) The proposed solvers based on the feed-forward neural network optimized with the GA,SQP,and GA-SQP algorithms can effectively solve different variants of nonlinear MHD Jeffery-Hamel problems once the original problem is transformed into an equivalently third-order boundary value problem of third-order ODEs.

(ii) The proposed solution of the SQP and GA-SQP methods matches with the standard solution up to 7 to 9 decimal places of accuracy, and is generally better than the reported results of the DTM,HPM, HAM,OHAM,and VIM for each case of both the problems.

(iii) In all the four cases of each Jeffery-Hamel problem,the results of the statistical analysis based on 100 independent runs for the proposed neural networks optimized with the GA,SQP,and GA-SQP techniques show that the values of the statistical parameters,i.e.,$G_{\rm {MIN}}$,$G_{\rm {MAE}}$,and $G_{\rm {STD}}$,fitness values,$G_{\rm MAE}$,and $M_{\rm F}$,are the best for the hybrid GA-SQP approach when they are compared with those of the GA and SQP algorithms. Generally,the most accurate and consistently convergent results are obtained with the GA-SQP method.

(iv) The comparative time analysis of the results shows that the GA-SQP algorithm takes comparatively longer time for the execution than that of the GA and SQP techniques. However,expected compensation can be given due to its invariable excellence in the accuracy and convergence observed for each case of both the MHD Jeffery-Hamel problems.

Appendix A

Here, a set of optimal weights are provided for the first four cases of both the Jeffery-Hamel problems. These weights are used in (15) to reproduce the proposed results of the designed approaches based on the GA, SQP, and GA-SQP algorithms.