*Article* **The Least Squares Homotopy Perturbation Method for Systems of Differential Equations with Application to a Blood Flow Model**

**M ˘ad ˘alina Sofia Pa¸sca 1,2,\*,†, Olivia Bund ˘au 1,†, Adina Juratoni 1,† and Bogdan C ˘aruntu 1,†**


† These authors contributed equally to this work.

**Abstract:** In this paper, least squares homotopy perturbation is presented as a straightforward and accurate method to compute approximate analytical solutions for systems of ordinary differential equations. The method is employed to solve a problem related to a laminar flow of a viscous fluid in a semi-porous channel, which may be used to model the blood flow through a blood vessel, taking into account the effects of a magnetic field. The numerical computations show that the method is both easy to use and very accurate compared to the other methods previously used to solve the given problem.

**Keywords:** least squares homotopy perturbation method; system of nonlinear differential equations; approximate analytical solutions; non-Newtonian fluid; magnetohydrodynamics

#### **1. Introduction**

The least squares homotopy perturbation method was introduced in 2017 by Bota and Caruntu in [1], and its main feature is an accelerated convergence compared to the regular homotopy perturbation method. In the few years since its introduction, the method (or slightly modified versions of it) was used by several researchers with very good results in finding approximate solutions for various types of problems, among which, we mention:


In the present paper, we employ the least squares homotopy perturbation method to compute approximate analytical solutions for boundary problems consisting of systems of nonlinear ordinary differential equations of the type:

$$\mathcal{L}\_i(\mathsf{U}\_1(y), \mathsf{U}\_2(y), \dots, \mathsf{U}\_n(y)) + \mathcal{N}\_i(\mathsf{U}\_1(y), \mathsf{U}\_2(y), \dots, \mathsf{U}\_n(y)) - f\_i(y) = 0, \qquad i = \overline{1, n} \tag{1}$$

$$\mathcal{B}\_j(\mathsf{U}\_i(y)) = 0, \quad j = \overline{1, k} \tag{1}$$

where *Ui*(*y*) are the unknown functions, L*<sup>i</sup>* are linear operators, N*<sup>i</sup>* are nonlinear operators, *f*(*t*)*<sup>i</sup>* are given functions, *y* denotes the variable, and B*<sup>j</sup>* are boundary operators.

#### **2. The Least Squares Homotopy Perturbation Method**

In this section, the least squares homotopy perturbation method (LSHPM) is presented. Since the numerical application considered in the following section only contains two equations, we introduce LSHPM for the case of a system consisting of two equations. Obviously, LSHPM can be easily generalized for systems consisting of as many equations as needed. We should also note that LSHPM works as well, if instead of B*j*(*Ui*(*y*)) = 0, we

**Citation:** Pa¸sca, M.S.; Bund ˘au, O.; Juratoni, A.; C ˘aruntu, B. The Least Squares Homotopy Perturbation Method for Systems of Differential Equations with Application to a Blood Flow Model. *Mathematics* **2022**, *10*, 546. https://doi.org/10.3390/ math10040546

Academic Editor: Marco Pedroni

Received: 2 November 2021 Accepted: 6 February 2022 Published: 10 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

have other types of conditions, such as B*j*(*Ui*(*y*), *Vi*(*y*)) = 0, initial-type conditions, or any combinations of the above.

Thus we consider the system:

$$\begin{array}{ll} \mathcal{L}\_1(\mathcal{U}(y), \mathcal{V}(y)) + \mathcal{N}\_1(\mathcal{U}(y), \mathcal{V}(y)) - f\_1(y) = 0, \\ \mathcal{L}\_2(\mathcal{U}(y), \mathcal{V}(y)) + \mathcal{N}\_2(\mathcal{U}(y), \mathcal{V}(y)) - f\_2(y) = 0, \\ \mathcal{B}\_i(\mathcal{U}(y)) = 0, \; \mathcal{B}\_j(\mathcal{V}(y)) = 0, \quad i = \overline{1, n1} \quad j = \overline{1, n2} \end{array} \tag{2}$$

where *U*(*y*) and *V*(*y*) are the unknown functions, L1,L<sup>1</sup> are linear operators, N1, N<sup>2</sup> are nonlinear operators, and B*i*,B*<sup>j</sup>* are boundary operators.

Let *U*˜ (*y*) and *V*˜(*y*) be approximate solutions of the system (2). The error obtained by replacing the exact solutions *U*(*y*), *V*(*y*) of the system (2) with the approximate ones *U*˜ (*y*); *V*˜(*y*) is given by the remainders:

$$\begin{cases} \mathcal{R}\_1(y, \hat{\mathcal{U}}) = \mathcal{L}\_1(\hat{\mathcal{U}}(y)) + \mathcal{N}\_1(\hat{\mathcal{U}}(y)) - f\_1(y), \ y \in \mathbb{R} \\ \mathcal{R}\_2(y, \hat{\mathcal{V}}) = \mathcal{L}\_2(\hat{\mathcal{V}}(y)) + \mathcal{N}\_2(\hat{\mathcal{V}}(y)) - f\_2(y), \ y \in \mathbb{R} \end{cases} \tag{3}$$

Following the homotopy perturbation method [9–11], the first step in applying LSHPM is to attach to the system (3) the family of equations:

$$\begin{bmatrix} (1-p)[\mathcal{L}\_1(\Phi\_1(y,p)) - f\_1(y)] + p\left[\mathcal{L}\_1(\Phi\_1(y,p)) + \mathcal{N}\_1(\Phi\_1(y,p)) - f\_1(y)\right] = 0\\ (1-p)[\mathcal{L}\_2(\Phi\_2(y,p)) - f\_2(y)] + p\left[\mathcal{L}\_2(\Phi\_2(y,p)) + \mathcal{N}\_2(\Phi\_2(y,p)) - f\_2(y)\right] = 0 \end{bmatrix} \tag{4}$$

where *p* ∈ [0, 1] is an embedding parameter and Φ*i*(*y*, *p*) with *i* = 1, 2 are unknown functions.

When *p* increases from 0 to 1, the solutions Φ*i*(*y*, *p*) of system (4) vary from Φ1(*y*, 0) = *U*0(*y*) and Φ2(*y*, 0) = *V*0(*y*) to the solutions Φ1(*y*, 1) = *U*(*y*) and Φ2(*y*, 1) = *V*(*y*) of the system (2). The functions *U*0(*y*) and *V*0(*y*) are the solutions of the system:

$$\begin{array}{c} \mathcal{L}\_1(\mathsf{U}\_0(y)) - f\_1(y) = 0 \\ \mathcal{L}\_2(V\_0(y)) - f\_2(y) = 0 \\ \mathcal{B}\_i(\mathsf{U}(y)) = 0, \; \mathcal{B}\_j(V(y)) = 0, \quad i = \overline{1, n1} \quad j = \overline{1, n2} \end{array} \tag{5}$$

We consider the following expansions of Φ*i*(*y*, *p*):

$$\begin{aligned} \Phi\_1(y, p) &= \mathcal{U}\_0(y) + \sum\_{m \ge 1} \mathcal{U}\_m(y) \ p^m \\ \Phi\_2(y, p) &= V\_0(y) + \sum\_{m \ge 1} V\_m(y) \ p^m \end{aligned} \tag{6}$$

Substituting the relations (6) in (4), collecting the same powers of *p* and equating the coefficients of the powers of *p*, we obtain:

$$\begin{array}{l} \mathcal{L}\_1(\mathcal{U}\_m(y)) = -\mathcal{N}\_1^{m-1}(\mathcal{U}\_0(y), \mathcal{U}\_1(y), \dots, \mathcal{U}\_{m-1}(y)) \\ \mathcal{L}\_2(\mathcal{V}\_m(y)) = -\mathcal{N}\_1^{m-1}(\mathcal{V}\_0(y), \mathcal{V}\_1(y), \dots, \mathcal{V}\_{m-1}(y)) \\ \mathcal{B}\_i(\mathcal{U}\_m(y)) = 0, \; \mathcal{B}\_j(\mathcal{V}\_m(y)) = 0, \quad i = \overline{1, n1} \quad j = \overline{1, n2} \end{array} \tag{7}$$

where <sup>N</sup> *<sup>j</sup> <sup>i</sup>* , *<sup>j</sup>* <sup>≥</sup> 0 are the coefficients of *<sup>p</sup><sup>j</sup>* in the nonlinear operator <sup>N</sup>*i*:

$$\begin{array}{c} \mathcal{N}\_{1}(\mathcal{U}(\boldsymbol{y})) = \mathcal{N}\_{1}^{0}(\mathcal{U}\_{0}(\boldsymbol{y})) + p\mathcal{N}\_{1}^{1}(\mathcal{U}\_{0}(\boldsymbol{y}), \mathcal{U}\_{1}(\boldsymbol{y})) + p^{2}\mathcal{N}\_{1}^{2}(\mathcal{U}\_{0}(\boldsymbol{y}), \mathcal{U}\_{1}(\boldsymbol{y}), \mathcal{U}\_{2}(\boldsymbol{y})) + \dots \\ \mathcal{N}\_{2}(\mathcal{V}(\boldsymbol{y})) = \mathcal{N}\_{1}^{0}(\mathcal{V}\_{0}(\boldsymbol{y})) + p\mathcal{N}\_{1}^{1}(\mathcal{V}\_{0}(\boldsymbol{y}), \mathcal{V}\_{1}(\boldsymbol{y})) + p^{2}\mathcal{N}\_{1}^{2}(\mathcal{V}\_{0}(\boldsymbol{y}), \mathcal{V}\_{1}(\boldsymbol{y}), \mathcal{V}\_{2}(\boldsymbol{y})) + \dots \end{array} \tag{8}$$

Now we can denote by

$$\begin{array}{l} f\_{1m} = \mathcal{U}\_0 + \mathcal{U}\_1 + \dots + \mathcal{U}\_{m\prime} \\ f\_{2m} = V\_0 + V\_1 + \dots + V\_m \end{array} \tag{9}$$

where *Um*, *m* ≥ 1, and *Vm*, *m* ≥ 1, are obtained from the linear Equation (7).

For *m* = 0, 1, 2, . . . let us consider the set *Sim* containing the functions

*ϕim*0, *ϕim*1, *ϕim*2, . . ., *ϕimnm* , (10)

chosen as linearly independent functions in the vector space of the continuous functions on the real interval *I*, such that *Sim*−<sup>1</sup> ⊆ *Sim* and *fim* is a real linear combination of these functions where *i* = 1, 2.

Using the functions given by (10), we define some types of approximate solutions of the system (2).

**Definition 1.** *A sequence of functions* {*sim*(*y*)}*m*∈<sup>N</sup> *of the form*

$$s\_{im}(\underline{y}) = \sum\_{k=0}^{n\_{m}} a\_{im}^{k} \varphi\_{imk} \; , \; m \in \mathbb{N} \; , \; a\_{m}^{k} \in \mathbb{R} \; ; \; i = \overline{1,2} \tag{11}$$

*are called HP-sequences of system (2).*

*Functions of the HP-sequences are called HP-functions of system (2). The HP-sequences* {*sim*(*y*)}*m*∈<sup>N</sup> *with the property*

$$\lim\_{m \to \infty} R\_i(y, \mathbf{s}\_{1m}(y), \mathbf{s}\_{2m}(y)) = 0, \quad i = \overline{1,2}$$

*is called convergent to the solution of the system (2).*

**Definition 2.** *The HP-functions U and* ˜ *V satisfying the conditions* ˜

<sup>|</sup>*R*1(*y*, *<sup>U</sup>*˜ , *<sup>V</sup>*˜)<sup>|</sup> <sup>&</sup>lt; , <sup>B</sup>*i*(*U*˜ ) = <sup>0</sup> <sup>|</sup>*R*2(*y*, *<sup>U</sup>*˜ , *<sup>V</sup>*˜)<sup>|</sup> <sup>&</sup>lt; , <sup>B</sup>*j*(*V*˜) = <sup>0</sup> (12)

*are called ε-approximate HP-solutions of the system (2).*

**Definition 3.** *HP-function U and* ˜ *V satisfying the conditions* ˜

$$\int\_{\mathbb{T}} \mathbb{R}^2\_1(y, \tilde{\mathcal{U}}, \tilde{\mathcal{V}}) dy \le \delta, \quad \mathcal{B}\_i(\tilde{\mathcal{y}}) = 0 \tag{13}$$

$$\int\_{\mathbb{T}} \mathcal{R}\_2^2(y, \vec{v}) dy \le \delta, \quad \mathcal{B}\_{\vec{\jmath}}(\vec{y}) = 0 \tag{14}$$

*are called weak δ-approximate HP-solutions of the system (2) on the real interval I.*

**Remark 1.** *It is easy to see that any ε-approximate HP-solution of the system (2) is also a weak approximate HP-solution. It follows that the set of weak approximate HP-solutions of the system (2) also contains the approximate HP-solutions of the system.*

The following theorem states the existence of weak approximate HP-solutions of the system (2) and furnishes the way to construct them.

**Theorem 1.** *The system (2) admits a sequence of weak approximate HP-solutions.*

**Proof.** The first step of the proof is to construct the HP-sequences {*sim*(*t*)}*m*∈N, *<sup>i</sup>* = 1, 2. Let us consider the approximate HP-solutions of the type:

$$\begin{aligned} \mathcal{U} &= \sum\_{k=0}^{n\_m} c\_m^k \rho\_{1mk\prime} \text{ where } m = 0, 1, \dots \text{ and } \\ \mathcal{V} &= \sum\_{k=0}^{n\_m} d\_m^k \rho\_{2mk\prime} \text{ where } m = 0, 1, \dots \end{aligned}$$

In the following, the unknown constants *c<sup>k</sup> <sup>m</sup>* and *d<sup>k</sup> <sup>m</sup> k* ∈ {0, 1, ... , *km*}, will be determined. Substituting the approximate solutions *U*˜ and *V*˜ in the system (2), one gets the expression:

$$\begin{aligned} \mathcal{R}\_1(y, d\_m^k) &:= R\_1(y, \check{U}) \\ \mathcal{R}\_2(y, d\_m^k) &:= R\_2(y, \check{V}). \end{aligned} \tag{15}$$

Attaching to the system (2) the following real functionals

$$\begin{aligned} J\_1(c\_m^k) &= \int \mathcal{R}\_1^2(y, c\_m^k) dy \\ J\_2(d\_m^k) &= \int \mathcal{R}\_2^2(y, d\_m^k) dy \end{aligned} \tag{16}$$

and imposing the boundary conditions, we can determine *l* ∈ N, *l* ≤ *m*, such that *cm* <sup>0</sup> , *<sup>c</sup><sup>m</sup>* <sup>1</sup> , ..., *<sup>c</sup><sup>m</sup> <sup>l</sup>* and *<sup>d</sup><sup>m</sup>* <sup>0</sup> , *<sup>d</sup><sup>m</sup>* <sup>1</sup> , ..., *<sup>d</sup><sup>m</sup> <sup>l</sup>* are computed as functions of *<sup>c</sup><sup>m</sup> <sup>l</sup>*+1, *<sup>c</sup><sup>m</sup> <sup>l</sup>*+2, ..., *<sup>c</sup><sup>m</sup> <sup>n</sup>* respectively *dm <sup>l</sup>*+1, *<sup>d</sup><sup>m</sup> <sup>l</sup>*+2, . . ., *<sup>d</sup><sup>m</sup> n* .

Replacing *c<sup>m</sup>* <sup>0</sup> , *<sup>c</sup><sup>m</sup>* <sup>1</sup> , ..., *<sup>c</sup><sup>m</sup> <sup>l</sup>* and *<sup>d</sup><sup>m</sup>* <sup>0</sup> , *<sup>d</sup><sup>m</sup>* <sup>1</sup> , ..., *<sup>d</sup><sup>m</sup> <sup>l</sup>* in (16), the values of *c*˜ *m <sup>l</sup>*+1, *c*˜ *m <sup>l</sup>*+2, ..., *c*˜ *m <sup>n</sup>* respectively ˜*d<sup>m</sup> <sup>l</sup>*+1, ˜*d<sup>m</sup> <sup>l</sup>*+2, ..., ˜*d<sup>m</sup> <sup>n</sup>* are computed as the values, which give the minimum of the functional (16).

Using again the boundary conditions, the values of *c*˜ *m* <sup>0</sup> , *c*˜ *m* <sup>1</sup> , . . ., *c*˜ *m <sup>l</sup>* as functions of *c*˜ *m <sup>l</sup>*+1, *c*˜ *m <sup>l</sup>*+2, ..., *c*˜ *m <sup>n</sup>* are determined and the values of ˜*d<sup>m</sup>* <sup>0</sup> , ˜*d<sup>m</sup>* <sup>1</sup> , ..., ˜*d<sup>m</sup> <sup>l</sup>* as functions of ˜*d<sup>m</sup> <sup>l</sup>*+1, ˜*d<sup>m</sup> <sup>l</sup>*+2, ..., ˜*d<sup>m</sup> n* are determined.

Using the constants *c*˜ *m* <sup>0</sup> , ..., *c*˜ *m <sup>n</sup>* and ˜*d<sup>m</sup>* <sup>0</sup> , ..., ˜*d<sup>m</sup> <sup>n</sup>* thus determined, the following HP-functions

$$\begin{aligned} s\_{1m}(t) &= \sum\_{k=0}^{n\_m} \vec{c}\_m^k \, \boldsymbol{\varrho}\_{mk} \\ s\_{2m}(t) &= \sum\_{k=0}^{n\_m} \vec{d}\_m^k \, \boldsymbol{\varrho}\_{mk} \end{aligned} \tag{17}$$

are constructed.

The second step of the proof is to show that the above HP-functions *sim*(*y*) are weak approximate solutions of the system (2).

Based on the way the HP-functions *sim*(*y*) are computed and taking into account that *fim* given by (9) are HP-functions for system (2), it follows:

$$0 \le \int\_I \mathcal{R}\_i^2(y, s\_{im}(y)) dy \le \int\_I \mathcal{R}\_i^2(y, f\_{im}(y)) dy \; \; \forall m \in \mathbb{N}, i = \overline{1, 2}.$$

Therefore,

$$0 \le \lim\_{m \to \infty} \int\_I R\_i^2(y, s\_{im}(y)) dy \le \lim\_{m \to \infty} \int\_I R\_i^2(y, f\_{im}(y)) dy, i = \overline{1,2}.$$

Since *fim* are convergent to the solution of the system (2), we obtain:

$$\lim\_{m \to \infty} \int\_I R\_i^2(y, s\_{im}(y)) dy = 0.$$

It follows that for all > 0 there exists *m*<sup>0</sup> ∈ N such that for all *m* ∈ N, *m* > *m*0, the sequence *sim*(*y*) is a weak -approximate HP-solution of the system (2).

**Remark 2.** *The proof of the above theorem give us a way to determine a weak approximate HPsolution of the system (2), <sup>U</sup>*˜ , *<sup>V</sup>*˜ *. Moreover, taking into account the Remark 1, if* <sup>|</sup>*R*1(*y*, *<sup>U</sup>*˜ )<sup>|</sup> <sup>&</sup>lt; *and* <sup>|</sup>*R*(*y*, *<sup>V</sup>*˜)<sup>|</sup> <sup>&</sup>lt; *then U and* ˜ *V are also* ˜ *-approximate HP-solutions of the considered system.*

#### **3. Numerical Application**

The application presented in this section is the one included in the paper by Basiri Parsa, Rashidi et al. [12], where the authors employed the well-known homotopy analysis method and differential transform method to find approximate analytical solutions for the following boundary value problem:

$$\begin{array}{c} \mathsf{U}V' - V\mathsf{U}' = \frac{1}{\mathsf{Re}} \left( \mathsf{U}' - Ha^2 \mathsf{U} \right) \\ V^{\mathsf{V}} = Ha^2 V'^\prime + \mathrm{Re}(V'V'^\prime - VV^{\prime\prime}) \\ \mathsf{U}(0) = 1, \; V(0) = 0, \; V'(0) = 0, \; \mathsf{U}(1) = 0, \; V(1) = 1, \; V'(1) = 0 \end{array} \tag{18}$$

These equations model a laminar magnetohydrodynamic flow of a non-Newtonian viscous fluid in a semi-porous channel under the influence of an axial uniform static magnetic field. *U* and *V* are the mean axial and normal velocity components, respectively, *Ha* is the Hartmann number, and *Re* is the Reynolds number.

In order to find analytical solutions for this type of problems, various approximation methods were employed over the years, with various rates of success, methods among which we mention: the homotopy perturbation method [10], the variational iteration method [11], the Adomian decomposition method [13], and the optimal homotopy asymptotic method [14]. While these methods (and many others) have been successfully employed, due to the nature of the equations, the computations involved are usually very difficult.

We remark that the system (18) may be used to study the influence of a magnetic field on the blood flow through a blood vessel.Numerous models have been established for the study of the hydrodynamic blood flow through the vessels, for example in [15]. Here, the authors analyze the flow of blood in tubes with reduced diameters, while in [16], the authors engage themselves in the study of the blood flow in small arched tubes, which are modeled. Blood flow has been analyzed trough the effect of the magnetic field as an great electrically conductive fluid. Knowing that blood is a ferrofluid, it can be concluded that there is the possibility of controlling the blood pressure and its flow behavior by using a fitting magnetic field. In [17], the authors came up with a mathematical representation of the blood flow in a blood vessel of reduced dimensions, in the presence of a magnetic field. Moreover, in [18], the authors investigated the apparatus of interaction between the red blood cells and an external magnetic field. The results will show the capacity of a magnetic field to modulate the blood flow. Other research on the magnetic properties of the blood are based on [12,19–27].

Many mathematical models show parts of the human circulatory system (for example [28–31]), most of the time, the blood flow is modeled by using differential equations, mostly nonlinear ones. However, it is usually nearly impossible to find exact solutions for these types of equations. Such cases require approximation methods for calculating almost exact solutions, because these approximated solutions may provide important information about the phenomenon.

In the following, we apply the least squares homotopy perturbation method to compute approximate polynomial solutions for the system (18) for two cases with particular significant values of the Hartmann number *Ha* and of the Reynolds number *Re*, and we compare our solutions with previous ones obtained in the literature.

#### *3.1. The Case Re* = 1 *and Ha* = 0

The case *Re* = 1 and *Ha* = 0 corresponds to a non-conducting blood flow. In [12], Basiri Parsa et al. computed approximate solutions of the system (18) by using the homotopy analysis method (HAM) and the differential transform method (DTM), and in [32], Caruntu et al. computed approximate solutions by using the polynomial least squares method (PLSM).

In this case, employing LSHPM for the system (18), we compute the approximate solutions as follows:

The linear operators are:

$$\begin{aligned} \mathcal{L}\_1(\Phi\_1(y, p)) &= -\frac{1}{R\mathfrak{e}\_1} \mathcal{U}^{\prime} \\ \mathcal{L}\_2(\Phi\_2(y, p)) &= \mathcal{V}^{\Gamma} \end{aligned} \tag{19}$$

and the nonlinear operators are:

$$\begin{split} \mathcal{N}\_1(\Phi\_1(t,p)) &= \frac{Ha^2}{\text{Re}} \mathcal{U} + \mathcal{U}V' - \mathcal{U}'V \\ \mathcal{N}\_2(\Phi\_2(t,p)) &= -Ha^2V' - \text{Re}(V'V'^\prime - VV''^\prime). \end{split} \tag{20}$$

We obtain the HPM approximations:

• First-term approximations:

$$\begin{array}{l} \mathcal{U}\_0(y) = 1 - y \\ \mathcal{V}\_0(y) = 3y^2 - y^3 \end{array} \tag{21}$$

• Second-term approximations:

$$\begin{array}{c} \mathcal{U}\_{1}(y) = \frac{y^{5}}{5} - \frac{3y^{4}}{4} + y^{3} - \frac{29y}{20} + 1\\ \mathcal{V}\_{1}(y) = \frac{2y^{7}}{35} - \frac{y^{6}}{5} + \frac{3y^{5}}{10} - \frac{167y^{3}}{70} + \frac{113y^{2}}{35} \end{array} \tag{22}$$

• Third-term approximations:

$$\begin{array}{c} \text{II}\_{2}(y) = \frac{2y^{9}}{315} - \frac{19y^{8}}{560} + \frac{y^{7}}{20} - \frac{y^{6}}{20} + \frac{23y^{5}}{70} - \frac{1643y^{4}}{1680} + \frac{113y^{3}}{105} - \frac{1763y}{1260} + 1\\ \text{Vr}(y) = \frac{4y^{11}}{5775} - \frac{2y^{10}}{525} + \frac{y^{9}}{210} - \frac{3y^{8}}{560} + \frac{97y^{7}}{1225} - \frac{533y^{6}}{2100} + \frac{121y^{5}}{350} - \frac{774469y^{3}}{323400} + \frac{2087479y^{2}}{646800} \end{array} \tag{23}$$

For the second-term approximation *U*<sup>0</sup> + *U*1, the set *S*1*<sup>m</sup>* consist of the functions {*y*, *<sup>y</sup>*3, *<sup>y</sup>*4, *<sup>y</sup>*5} and for *<sup>V</sup>*<sup>0</sup> <sup>+</sup> *<sup>V</sup>*<sup>1</sup> the set is *<sup>S</sup>*2*<sup>m</sup>* <sup>=</sup> {*y*2, *<sup>y</sup>*3, *<sup>y</sup>*5, *<sup>y</sup>*6, *<sup>y</sup>*7}.

We will compute the approximate solutions *U*˜ (*y*) = *c*<sup>0</sup> + *c*1*y* + *c*2*y*<sup>3</sup> + *c*3*y*<sup>4</sup> + *c*4*y*<sup>5</sup> and *V*˜(*y*) = *d*0*y*<sup>2</sup> + *d*1*y*<sup>3</sup> + *d*2*y*<sup>5</sup> + *d*3*y*<sup>6</sup> + *d*4*y*7.

From the initial conditions: *U*˜ (0) = 1, *U*˜ (1) = 0 and *V*˜(0) = 0, *V*˜(1) = 1, *V*˜  (0) = 0, *V*˜  (1) = 0 we obtain: *c*<sup>0</sup> = 1 and *c*<sup>1</sup> = −1 − *c*<sup>2</sup> − *c*<sup>3</sup> − *c*<sup>4</sup> respectively −*d*<sup>0</sup> = 2*d*<sup>2</sup> + 3*d*<sup>3</sup> + 4*d*<sup>4</sup> + 3 and *d*<sup>1</sup> = −3*d*<sup>2</sup> − 4*d*<sup>3</sup> − 5*d*<sup>4</sup> − 2.

Replacing these expressions of *c*0, *c*1, *d*0, and *d*<sup>1</sup> in the corresponding remainders (15) are:

$$\begin{aligned} \mathcal{R}\_1(y,\Omega) &= \mathcal{R}(y, c\_2, c\_3, c\_4) \\ \mathcal{R}\_2(y,\bar{V}) &= \mathcal{R}(y, d\_2, d\_3, d\_4) \end{aligned} \tag{24}$$

Next, we compute the corresponding functionals (16) (too long to be included here):

$$\begin{aligned} f\_1(c\_2, c\_3, c\_4) &= \int \mathcal{R}\_1^2(\mathcal{y}, c\_2, c\_3, c\_4) d\mathcal{y} \\ f\_2(d\_2, d\_3, d\_4) &= \int \mathcal{R}\_2^2(\mathcal{y}, d\_2, d\_3, d\_4) d\mathcal{y} \end{aligned} \tag{25}$$

and we compute the minimum of this functionals, determining the coefficients *cj* and *dj*, *j* = 2, 4, thus finding the approximate solutions of the system (18) by LSHPM.

In a similar way, we compute the third-term approximations by LSHPM, and the solutions are:

	- *<sup>U</sup>*˜ (*y*) = <sup>−</sup>0.012783174216369037776*y*<sup>9</sup> <sup>+</sup> 0.079312017411340841618*y*<sup>8</sup> <sup>−</sup> 0.21264602372589301173*y*<sup>7</sup> <sup>+</sup> 0.26096763633169609359*y*<sup>6</sup>

<sup>+</sup> 0.13236850157788926145*y*<sup>5</sup> <sup>−</sup> 0.90665763559833919190*y*<sup>4</sup> <sup>+</sup> 1.0653323749944480412*y*<sup>3</sup> − 1.4058936967747729964*y* + 1


The comparison of these LSHPM solutions with the previous approximate solutions computed in [12] using HAM and DTM, and in [32] using PLSM, is presented in Figures 1 and 2. Since no exact solutions are available, the comparison is done by means of computing the error relative to a fourth-order Runge–Kutta method numerical solution (i.e., the absolute errors are computed as the difference between our approximate solutions and the numerical solutions).

**Figure 1.** Comparison of absolute errors corresponding to the approximation from [12] *UDTM* (red curve), [10] *UHPM* 3 terms (blue curve), [32] *UPLSM* (orange curve), and our LSHPM approximation *ULSHPM* (green curve).

**Figure 2.** Comparison of absolute errors corresponding to the approximation from [12] *VDTM* (red curve), [10] *VHPM* 3 terms (blue curve), [32] *VPLSM* (orange curve), and our LSHPM approximation *VLSHPM* (green curve).

The comparison is further illustrated by Tables 1 and 2, which includes the results obtained in [12], by means of the HAM and DTM, the results obtained in [32] by PLSM, and the results computed by classical HPM and by LSHPM. The comparison lead to the same conclusion: the approximate solutions obtained by LSHPM are more precise.


**Table 1.** Comparison of the absolute errors of the approximate solutions *U* in case *Re* = 1 and *Ha* = 0.

**Table 2.** Comparison of the absolute errors of the approximate solutions *V* in case *Re* = 1 and *Ha* = 0.


*3.2. The Case Re* = 1 *and Ha* = 1

In the case *Re* = 1 and *Ha* = 1, the influence of the magnetic field on the blood flow is non-negligible and the flow is weakly magnetic. The computations by LSHPM are similar to the ones in the previous case.

The approximations terms by HPM are:

• First-term approximations:

$$\begin{array}{l} \mathcal{U}\_0(y) = 1 - y \\ \mathcal{V}\_0(y) = 3y^2 - 2y^3 \end{array} \tag{26}$$

• Second-term approximations:

$$\begin{array}{c} \mathcal{U}\_{1}(y) = \frac{y^{5}}{5} - \frac{3y^{4}}{4} + \frac{5y^{3}}{6} + \frac{y^{2}}{2} - \frac{107y}{60} + 1\\\mathcal{V}\_{1}(y) = \frac{2y^{7}}{35} - \frac{y^{6}}{5} + \frac{y^{5}}{5} + \frac{y^{4}}{4} - \frac{181y^{3}}{70} + \frac{459y^{2}}{140} \end{array} \tag{27}$$

• Third-term approximations:

$$\begin{array}{c} \text{U2}(y) = \frac{2y^9}{315} - \frac{19y^8}{560} + \frac{9y^7}{140} - \frac{2y^6}{15} + \frac{2129y^5}{4200} - \frac{451y^4}{420} + \frac{401y^3}{504} + \frac{y^2}{2} - \frac{41129y}{25200} + 1\\\text{V2}(y) = \frac{4y^{11}}{5775} - \frac{2y^{10}}{525} + \frac{y^9}{140} - \frac{9y^8}{560} + \frac{1507y^7}{14700} - \frac{1129y^6}{4200} + \frac{317y^5}{1400} + \frac{153y^4}{560}\\\ \phantom{\text{V2}} - \frac{838379y^3}{323400} + \frac{352823y^2}{107800} \end{array} \tag{28}$$

The corresponding solutions obtained by using LSHPM are:

• Second-term approximations:

*<sup>U</sup>*˜ (*y*) = 0.24252310472551437387*y*<sup>5</sup> <sup>−</sup> 0.79607170987751279666*y*<sup>4</sup>

	- *<sup>U</sup>*˜ (*y*) = <sup>−</sup>0.039570918805227895983*y*<sup>9</sup> <sup>+</sup> 0.18039801278862870509*y*<sup>8</sup>
	- <sup>−</sup>0.33738809541609613996*y*<sup>7</sup> <sup>+</sup>0.25556139173308396128*y*<sup>6</sup> <sup>+</sup>0.28898499486097221037*y*<sup>5</sup> <sup>−</sup>0.98687905924276726841*y*<sup>4</sup> <sup>+</sup>0.80192640217575596712*y*<sup>3</sup> <sup>+</sup>0.50063075107844706510*y*<sup>2</sup>
	- − 1.6636634791727966046*y* + 1
	- *<sup>V</sup>*˜(*y*) = 0.00070249628967946626576*y*<sup>11</sup> <sup>−</sup> 0.0034005002946346048030*y*<sup>10</sup>
	- <sup>+</sup> 0.0063398061612668565765*y*<sup>9</sup> <sup>−</sup> 0.017261143110682215908*y*<sup>8</sup>
	- <sup>+</sup>0.10727677845680360313*y*<sup>7</sup> <sup>−</sup>0.27266199870132195652*y*<sup>6</sup> <sup>+</sup>0.22694957760601929024*y*<sup>5</sup>
	- <sup>+</sup> 0.27260900187074289280*y*<sup>4</sup> <sup>−</sup> 2.5917328827530869045*y*<sup>3</sup> <sup>+</sup> 3.2711788644752135727*y*<sup>2</sup>

Again, the comparisons included in Tables 3 and 4 allow us the reach the same conclusion as in the previous case, namely that the approximations obtained by LSHPM are more accurate than the previous ones by other methods.

We mention the fact that we computed approximations for a wide range of values of the parameters *Re* and *Ha*, and the above conclusions remained valid for all of the computed solutions.

**Table 3.** Comparison of the absolute errors of the approximate solutions for *U* for the case *Re* = 1 and *Ha* = 1.


**Table 4.** Comparison of the absolute errors of the approximate solutions for *V* for the case *Re* = 1 and *Ha* = 1.


#### **4. Discussion of the Results**

As we mention in the previous section, we computed approximations for a wide range of values of the parameters *Re* and *Ha*, and the conclusions of our study are in very good agreement with previous ones included in [12,33–35] and other studies. This is to be expected because, even though our approximations are much more precise than the ones included in the previous studies, the overall phenomena described by the solutions are obviously the same. In the following, we will briefly summarize the results of the study.

The first results are related to the influence of the Reynolds number on the velocities *U* and *V*, when the value of the Hartmann number is fixed. We were able to draw similar conclusions for both cases studied in [12] (the non-conducting case *Ha* = 0 and the weakly magnetic flow case *Ha* = 1). The consequences of any increase of the Reynolds number are a modest increase to the *V*(*y*) component of the velocity of the blood flow, and a major decrease of the *U*(*y*) component. The effect of this phenomenon is a major deceleration of the blood velocity in the x-direction. For the case *Ha* = 0, these conclusions are illustrated by the Figures 3 and 4, while for the case *Ha* = 1 (and, actually, for any other value of *Ha* on its nominal interval [0, 2]), the figures look very similar.

**Figure 3.** The effect of the increase of *Re* on *U* for the case *Ha* = 0—three dimensional plot.

**Figure 4.** The effect of the increase of *Re* on *V* for the case *Ha* = 0—three dimensional plot.

The next study item is the influence of the Hartmann number on the velocities *U* and *V* for fixed values of the Reynolds number. Because the magnetic field is applied in the y-direction only, there is no visible influence of the increase of the *Ha V*(*y*) component of the blood velocity while there is a decrease of the *U*(*y*) component, as expected. These conclusions are illustrated for *Re* = 1 in Figures 5 and 6 (again, we mention that for any values of the Reynolds number on its nominal interval [1, 20], the figures are similar).

**Figure 5.** The effect of the increase of *Ha* on *U* for the case *Re* = 1—three dimensional plot.

**Figure 6.** The effect of the increase of *Ha* on *V* for the case *Re* = 1—three dimensional plot.

In the last part of the study, we investigated the merged impact of *Re* and *Ha* on *U*(*y*) and *V*(*y*), impact highlighted in the Figures 7 and 8 for *y* = 0.2 (red surface), *y* = 0.4 (blue surface), *y* = 0.6 (yellow surface) and *y* = 0.8 (green surface).

Figure 7 is a good synthesis of the research done on the impact of *Re* and *Ha* on a replicated blood flow in a semi-porous channel, as it is obvious that the increase in both *Re* and *Ha* conduct to the reduction of *U*(*y*). At the same time, Figure 8 gives new insight regarding the interplay between the impact of *Re* and *Ha* on the blood flow velocity in the y-direction, particularly the fact that the decelerative consequence of the increase of *Ha* greatly depends on *Re*. This effect is greater for reduced values of *Re* and, as a result, in situations where, due to practical discussions, a strong suction at the upper wall ( defined by a large value of *Re*) cannot be achieved, a decrease of blood flow velocity can be achieved by boosting the intensity of the magnetic field applied. Even if it is practically possible, an increase of the suction at the upper wall is apparently the preferable method for reducing the flow, since the effect of the increase of Re seems to be considerably greater than the effect of the increase of Ha. Furthermore, if the value of the Reynolds number is large, the consequences of the increase of the Hartmann number is small.

**Figure 7.** The combined influence of *Re* and *Ha* on *U*. The red surface corresponds to *y* = 0.2, the blue surface to *y* = 0.4, the yellow surface to *y* = 0.6 and the green one to *y* = 0.8.

**Figure 8.** The combined influence of *Re* and *Ha* on *V*. The red surface corresponds to *y* = 0.2, the blue surface to *y* = 0.4, the yellow surface to *y* = 0.6, and the green one to *y* = 0.8.

#### **5. Conclusions**

The least squares homotopy perturbation method is introduced as a straightforward and very accurate method to compute analytically approximate solutions for systems of ordinary differential equations.

The method was employed for a system of equations modeling the blood flow through a blood vessel under the action of a magnetic field. The comparison with approximate solutions computed by using well-known methods, such as the homotopy analysis method, the differential transform method and the homotopy perturbation method, clearly illustrate the precision of our method.

**Author Contributions:** Conceptualization, M.S.P. and B.C.; Data curation, M.S.P., O.B., A.J. and B.C.; Formal analysis, B.C.; Funding acquisition, M.S.P., O.B., A.J. and B.C.; Investigation, M.S.P., O.B., A.J. and B.C.; Methodology, B.C.; Software, M.S.P. and B.C.; Supervision, B.C.; Validation, M.S.P. and B.C.; Visualization, M.S.P., A.J. and B.C.; Writing—original draft, M.S.P., O.B., A.J. and B.C.; Writing—review and editing, M.S.P. and B.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

