Next Article in Journal
A Linear, Second-Order, and Unconditionally Energy-Stable Method for the L2-Gradient Flow-Based Phase-Field Crystal Equation
Next Article in Special Issue
Statistical Analysis of Current Financial Instrument Quotes in the Conditions of Market Chaos
Previous Article in Journal / Special Issue
Lie Geometric Methods in the Study of Driftless Control Affine Systems with Holonomic Distribution and Economic Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Mădălina Sofia Paşca
1,2,*,†,
Olivia Bundău
1,†,
Adina Juratoni
1,† and
Bogdan Căruntu
1,†
1
Department of Mathematics, Politehnica University Timişoara, 300006 Timişoara, Romania
2
Department of Mathematics, West University of Timişoara, 300223 Timişoara, Romania
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2022, 10(4), 546; https://doi.org/10.3390/math10040546
Submission received: 2 November 2021 / Revised: 30 January 2022 / Accepted: 6 February 2022 / Published: 10 February 2022

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.

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:
  • Boundary value problems for ordinary differential equations [2,3].
  • Fractional partial differential equations [4,5,6].
  • Fractional order integro-differential equations [7].
  • Systems of fractional partial differential equations [8].
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:
L i ( U 1 ( y ) , U 2 ( y ) , , U n ( y ) ) + N i ( U 1 ( y ) , U 2 ( y ) , , U n ( y ) ) f i ( y ) = 0 , i = 1 , n ¯ B j ( U i ( y ) ) = 0 , j = 1 , k ¯
where U i ( y ) are the unknown functions, L i are linear operators, N i are nonlinear operators, f ( t ) i are given functions, y denotes the variable, and B j 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 ( U i ( y ) ) = 0 , we have other types of conditions, such as B j ( U i ( y ) , V i ( y ) ) = 0 , initial-type conditions, or any combinations of the above.
Thus we consider the system:
L 1 ( U ( y ) , V ( y ) ) + N 1 ( U ( y ) , V ( y ) ) f 1 ( y ) = 0 , L 2 ( U ( y ) , V ( y ) ) + N 2 ( U ( y ) , V ( y ) ) f 2 ( y ) = 0 , B i ( U ( y ) ) = 0 , B j ( V ( y ) ) = 0 , i = 1 , n 1 ¯ j = 1 , n 2 ¯
where U ( y ) and V ( y ) are the unknown functions, L 1 , L 1 are linear operators, N 1 , N 2 are nonlinear operators, and B i , B j 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:
R 1 ( y , U ˜ ) = L 1 ( U ˜ ( y ) ) + N 1 ( U ˜ ( y ) ) f 1 ( y ) , y R R 2 ( y , V ˜ ) = L 2 ( V ˜ ( y ) ) + N 2 ( V ˜ ( y ) ) f 2 ( y ) , y R
Following the homotopy perturbation method [9,10,11], the first step in applying LSHPM is to attach to the system (3) the family of equations:
( 1 p ) [ L 1 ( Φ 1 ( y , p ) ) f 1 ( y ) ] + p [ L 1 ( Φ 1 ( y , p ) ) + N 1 ( Φ 1 ( y , p ) ) f 1 ( y ) ] = 0 ( 1 p ) [ L 2 ( Φ 2 ( y , p ) ) f 2 ( y ) ] + p [ L 2 ( Φ 2 ( y , p ) ) + N 2 ( Φ 2 ( y , p ) ) f 2 ( y ) ] = 0
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:
L 1 ( U 0 ( y ) ) f 1 ( y ) = 0 L 2 ( V 0 ( y ) ) f 2 ( y ) = 0 B i ( U ( y ) ) = 0 , B j ( V ( y ) ) = 0 , i = 1 , n 1 ¯ j = 1 , n 2 ¯
We consider the following expansions of Φ i ( y , p ) :
Φ 1 ( y , p ) = U 0 ( y ) + m 1 U m ( y ) p m Φ 2 ( y , p ) = V 0 ( y ) + m 1 V m ( y ) p m
Substituting the relations (6) in (4), collecting the same powers of p and equating the coefficients of the powers of p, we obtain:
L 1 ( U m ( y ) ) = N 1 m 1 ( U 0 ( y ) , U 1 ( y ) , , U m 1 ( y ) ) L 2 ( V m ( y ) ) = N 1 m 1 ( V 0 ( y ) , V 1 ( y ) , , V m 1 ( y ) ) B i ( U m ( y ) ) = 0 , B j ( V m ( y ) ) = 0 , i = 1 , n 1 ¯ j = 1 , n 2 ¯
where N i j , j 0 are the coefficients of p j in the nonlinear operator N i :
N 1 ( U ( y ) ) = N 1 0 ( U 0 ( y ) ) + p N 1 1 ( U 0 ( y ) , U 1 ( y ) ) + p 2 N 1 2 ( U 0 ( y ) , U 1 ( y ) , U 2 ( y ) ) + N 2 ( V ( y ) ) = N 1 0 ( V 0 ( y ) ) + p N 1 1 ( V 0 ( y ) , V 1 ( y ) ) + p 2 N 1 2 ( V 0 ( y ) , V 1 ( y ) , V 2 ( y ) ) +
Now we can denote by
f 1 m = U 0 + U 1 + + U m , f 2 m = V 0 + V 1 + + V m
where U m , m 1 , and V m , m 1 , are obtained from the linear Equation (7).
For m = 0 , 1 , 2 , let us consider the set S i m containing the functions
φ i m 0 , φ i m 1 , φ i m 2 , , φ i m n m ,
chosen as linearly independent functions in the vector space of the continuous functions on the real interval I, such that S i m 1 S i m and f i m 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 { s i m ( y ) } m N of the form
s i m ( y ) = k = 0 n m α i m k φ i m k , m N , α m k R , i = 1 , 2 ¯
are called HP-sequences of system (2).
Functions of the HP-sequences are called HP-functions of system (2).
The HP-sequences { s i m ( y ) } m N with the property
lim m R i ( y , s 1 m ( y ) , s 2 m ( y ) ) = 0 , i = 1 , 2 ¯
is called convergent to the solution of the system (2).
Definition 2.
The HP-functions U ˜ and V ˜ satisfying the conditions
| R 1 ( y , U ˜ , V ˜ ) | < ϵ , B i ( U ˜ ) = 0 | R 2 ( y , U ˜ , V ˜ ) | < ϵ , B j ( V ˜ ) = 0
are called ε-approximate HP-solutions of the system (2).
Definition 3.
HP-function U ˜ and V ˜ satisfying the conditions
I R 1 2 ( y , U ˜ , V ˜ ) dy δ , B i ( y ˜ ) = 0
I R 2 2 ( y , v ˜ ) dy δ , B j ( y ˜ ) = 0
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 { s i m ( t ) } m N , i = 1 , 2 ¯ .
Let us consider the approximate HP-solutions of the type:
U ˜ = k = 0 n m c m k φ 1 m k , where m = 0 , 1 , and
V ˜ = k = 0 n m d m k φ 2 m k , where m = 0 , 1 , .
In the following, the unknown constants c m k and d m k k { 0 , 1 , , k m } , will be determined.
Substituting the approximate solutions U ˜ and V ˜ in the system (2), one gets the expression:
R 1 ( y , c m k ) : = R 1 ( y , U ˜ ) R 2 ( y , d m k ) : = R 2 ( y , V ˜ ) .
Attaching to the system (2) the following real functionals
J 1 ( c m k ) = I R 1 2 ( y , c m k ) dy J 2 ( d m k ) = I R 2 2 ( y , d m k ) dy
and imposing the boundary conditions, we can determine l N , l m , such that c 0 m , c 1 m , , c l m and d 0 m , d 1 m , , d l m are computed as functions of c l + 1 m , c l + 2 m , , c n m respectively d l + 1 m , d l + 2 m , , d n m .
Replacing c 0 m , c 1 m , , c l m and d 0 m , d 1 m , , d l m in (16), the values of c ˜ l + 1 m , c ˜ l + 2 m , , c ˜ n m respectively d ˜ l + 1 m , d ˜ l + 2 m , , d ˜ n m are computed as the values, which give the minimum of the functional (16).
Using again the boundary conditions, the values of c ˜ 0 m , c ˜ 1 m , , c ˜ l m as functions of c ˜ l + 1 m , c ˜ l + 2 m , , c ˜ n m are determined and the values of d ˜ 0 m , d ˜ 1 m , , d ˜ l m as functions of d ˜ l + 1 m , d ˜ l + 2 m , , d ˜ n m are determined.
Using the constants c ˜ 0 m , , c ˜ n m and d ˜ 0 m , , d ˜ n m thus determined, the following HP-functions
s 1 m ( t ) = k = 0 n m c ˜ m k φ m k s 2 m ( t ) = k = 0 n m d ˜ m k φ m k
are constructed.
The second step of the proof is to show that the above HP-functions s i m ( y ) are weak approximate solutions of the system (2).
Based on the way the HP-functions s i m ( y ) are computed and taking into account that f i m given by (9) are HP-functions for system (2), it follows:
0 I R i 2 ( y , s i m ( y ) ) dy I R i 2 ( y , f i m ( y ) ) dy , m N , i = 1 , 2 ¯ .
Therefore,
0 lim m I R i 2 ( y , s i m ( y ) ) dy lim m I R i 2 ( y , f i m ( y ) ) dy , i = 1 , 2 ¯ .
Since f i m are convergent to the solution of the system (2), we obtain:
lim m I R i 2 ( y , s i m ( y ) ) dy = 0 .
It follows that for all ϵ > 0 there exists m 0 N such that for all m N , m > m 0 , the sequence s i m ( 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 HP-solution of the system (2), U ˜ , V ˜ . Moreover, taking into account the Remark 1, if | R 1 ( y , U ˜ ) | < ϵ and | R ( y , V ˜ ) | < ϵ 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:
U V V U = 1 R e U H a 2 U V I V = H a 2 V + R e V V V V U ( 0 ) = 1 , V ( 0 ) = 0 , V ( 0 ) = 0 , U ( 1 ) = 0 , V ( 1 ) = 1 , V ( 1 ) = 0
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, H a is the Hartmann number, and R e 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,20,21,22,23,24,25,26,27].
Many mathematical models show parts of the human circulatory system (for example [28,29,30,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 H a and of the Reynolds number R e , and we compare our solutions with previous ones obtained in the literature.

3.1. The Case R e = 1 and H a = 0

The case R e = 1 and H a = 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:
L 1 ( Φ 1 ( y , p ) ) = 1 R e U L 2 ( Φ 2 ( y , p ) ) = V I V
and the nonlinear operators are:
N 1 ( Φ 1 ( t , p ) ) = H a 2 R e U + U V U V N 2 ( Φ 2 ( t , p ) ) = H a 2 V R e V V V V .
We obtain the HPM approximations:
  • First-term approximations:
    U 0 ( y ) = 1 y V 0 ( y ) = 3 y 2 y 3
  • Second-term approximations:
    U 1 ( y ) = y 5 5 3 y 4 4 + y 3 29 y 20 + 1 V 1 ( y ) = 2 y 7 35 y 6 5 + 3 y 5 10 167 y 3 70 + 113 y 2 35
  • Third-term approximations:
    U 2 ( y ) = 2 y 9 315 19 y 8 560 + y 7 20 y 6 20 + 23 y 5 70 1643 y 4 1680 + 113 y 3 105 1763 y 1260 + 1 V 2 ( y ) = 4 y 11 5775 2 y 10 525 + y 9 210 3 y 8 560 + 97 y 7 1225 533 y 6 2100 + 121 y 5 350 774469 y 3 323400 + 2087479 y 2 646800
For the second-term approximation U 0 + U 1 , the set S 1 m consist of the functions { y , y 3 , y 4 , y 5 } and for V 0 + V 1 the set is S 2 m = { y 2 , y 3 , y 5 , y 6 , y 7 } .
We will compute the approximate solutions U ˜ ( y ) = c 0 + c 1 y + c 2 y 3 + c 3 y 4 + c 4 y 5 and V ˜ ( y ) = d 0 y 2 + d 1 y 3 + d 2 y 5 + d 3 y 6 + 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 0 = 1 and c 1 = 1 c 2 c 3 c 4 respectively d 0 = 2 d 2 + 3 d 3 + 4 d 4 + 3 and d 1 = 3 d 2 4 d 3 5 d 4 2 .
Replacing these expressions of c 0 , c 1 , d 0 , and d 1 in the corresponding remainders (15) are:
R 1 ( y , U ˜ ) = R ( y , c 2 , c 3 , c 4 ) R 2 ( y , V ˜ ) = R ( y , d 2 , d 3 , d 4 )
Next, we compute the corresponding functionals (16) (too long to be included here):
J 1 ( c 2 , c 3 , c 4 ) = 0 1 R 1 2 ( y , c 2 , c 3 , c 4 ) dy J 2 ( d 2 , d 3 , d 4 ) = 0 1 R 2 2 ( y , d 2 , d 3 , d 4 ) dy
and we compute the minimum of this functionals, determining the coefficients c j and d j , 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:
  • Second-term approximations:
    U ˜ ( y ) = 0.25139431098009168490 y 5 0.89052033909540654567 y 4
    + 1.0440071928069991650 y 3 1.4048811646916843042 y + 1
    V ˜ ( y ) = 0.058317127462336779863 y 7 0.21889392791262717538 y 6
    + 0.32522212582352969927 y 5 2.3916763031317642956 y 3 + 3.2270309777585249919 y 2
  • Third-term approximations:
    U ˜ ( y ) = 0.012783174216369037776 y 9 + 0.079312017411340841618 y 8
    0.21264602372589301173 y 7 + 0.26096763633169609359 y 6
    + 0.13236850157788926145 y 5 0.90665763559833919190 y 4 + 1.0653323749944480412 y 3
    1.4058936967747729964 y + 1
    V ˜ ( y ) = 0.00077252197269867863129 y 11 0.0033512835587865196224 y 10
    + 0.0028365610269174046158 y 9 0.0031087019658224288029 y 8
    + 0.079915383004597381009 y 7 0.25695066296630124619 y 6
    + 0.34706280997410686745 y 5 2.3943088377575857329 y 3 + 3.2271322102701755958 y 2
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 Figure 1 and Figure 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).
The comparison is further illustrated by Table 1 and Table 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.

3.2. The Case R e = 1 and H a = 1

In the case R e = 1 and H a = 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:
    U 0 ( y ) = 1 y V 0 ( y ) = 3 y 2 2 y 3
  • Second-term approximations:
    U 1 ( y ) = y 5 5 3 y 4 4 + 5 y 3 6 + y 2 2 107 y 60 + 1 V 1 ( y ) = 2 y 7 35 y 6 5 + y 5 5 + y 4 4 181 y 3 70 + 459 y 2 140
  • Third-term approximations:
    U 2 ( y ) = 2 y 9 315 19 y 8 560 + 9 y 7 140 2 y 6 15 + 2129 y 5 4200 451 y 4 420 + 401 y 3 504 + y 2 2 41129 y 25200 + 1 V 2 ( y ) = 4 y 11 5775 2 y 10 525 + y 9 140 9 y 8 560 + 1507 y 7 14700 1129 y 6 4200 + 317 y 5 1400 + 153 y 4 560 838379 y 3 323400 + 352623 y 2 107800
The corresponding solutions obtained by using LSHPM are:
  • Second-term approximations:
    U ˜ ( y ) = 0.24252310472551437387 y 5 0.79607170987751279666 y 4
    + 0.70039529434560076978 y 3 + 0.51680281524640722276 y 2
    1.6636495044400095698 y + 1
    V ˜ ( y ) = 0.058199385875154810585 y 7 0.20176414373359385855 y 6
    + 0.17947791112252870094 y 5 + 0.28462933749944678084 y 4 2.5916327628078782832 y 3
    + 3.2710902720443418494 y 2
  • Third-term approximations:
    U ˜ ( y ) = 0.039570918805227895983 y 9 + 0.18039801278862870509 y 8
    0.33738809541609613996 y 7 + 0.25556139173308396128 y 6 + 0.28898499486097221037 y 5
    0.98687905924276726841 y 4 + 0.80192640217575596712 y 3 + 0.50063075107844706510 y 2
    1.6636634791727966046 y + 1
    V ˜ ( y ) = 0.00070249628967946626576 y 11 0.0034005002946346048030 y 10
    + 0.0063398061612668565765 y 9 0.017261143110682215908 y 8
    + 0.10727677845680360313 y 7 0.27266199870132195652 y 6 + 0.22694957760601929024 y 5
    + 0.27260900187074289280 y 4 2.5917328827530869045 y 3 + 3.2711788644752135727 y 2
Again, the comparisons included in Table 3 and Table 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 R e and H a , and the above conclusions remained valid for all of the computed solutions.

4. Discussion of the Results

As we mention in the previous section, we computed approximations for a wide range of values of the parameters R e and H a , and the conclusions of our study are in very good agreement with previous ones included in [12,33,34,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 H a = 0 and the weakly magnetic flow case H a = 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 H a = 0 , these conclusions are illustrated by the Figure 3 and Figure 4, while for the case H a = 1 (and, actually, for any other value of H a on its nominal interval [0, 2]), the figures look very similar.
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 H a V ( y ) component of the blood velocity while there is a decrease of the U ( y ) component, as expected. These conclusions are illustrated for R e = 1 in Figure 5 and Figure 6 (again, we mention that for any values of the Reynolds number on its nominal interval [1, 20], the figures are similar).
In the last part of the study, we investigated the merged impact of R e and H a on U ( y ) and V ( y ) , impact highlighted in the Figure 7 and Figure 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 R e and H a on a replicated blood flow in a semi-porous channel, as it is obvious that the increase in both R e and H a conduct to the reduction of U ( y ) . At the same time, Figure 8 gives new insight regarding the interplay between the impact of R e and H a on the blood flow velocity in the y-direction, particularly the fact that the decelerative consequence of the increase of H a greatly depends on R e . This effect is greater for reduced values of R e and, as a result, in situations where, due to practical discussions, a strong suction at the upper wall ( defined by a large value of R e ) 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.

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

  1. Bota, C.; Caruntu, B. Approximate analytical solutions of nonlinear differential equations using the Least Squares Homotopy Perturbation Method. J. Math. Anal. Appl. 2017, 448, 401–408. [Google Scholar] [CrossRef]
  2. Bota, C.; Caruntu, B.; Lazureanu, C. The Least Squares Homotopy Perturbation Method for boundary value problems. Appl. Comput. Math. 2017, 16, 39–47. [Google Scholar]
  3. Qayyum, M.; Oscar, I. Least Square Homotopy Perturbation Method for Ordinary Differential Equations. J. Math. 2021, 2021, 7059194. [Google Scholar] [CrossRef]
  4. Thabet, H.; Kendre, S. Modified least squares homotopy perturbation method for solving fractional partial differential equations. Malaya J. Mat. 2018, 6, 420–427. [Google Scholar] [CrossRef] [Green Version]
  5. Kumar, R.; Koundal, R.; Shehzad, S.A. Generalized least square homotopy perturbation solution of fractional telegraph equations. Comput. Appl. Math. 2019, 38, 184. [Google Scholar] [CrossRef]
  6. Zhang, J.; Wei, Z.; Li, L.; Zhou, C. Least-Squares Residual Power Series Method for the Time-Fractional Differential Equations. Complexity 2019, 2019, 6159024. [Google Scholar] [CrossRef]
  7. Das, P.; Rana, S. Theoretical prospects of fractional order weakly singular Volterra Integro differential equations and their approximations with convergence analysis. Math. Methods Appl. Sci. 2021, 44, 9419–9440. [Google Scholar] [CrossRef]
  8. Kumar, R.; Koundal, R.; Shehzad, S.A. Modified homotopy perturbation approach for the system of fractional partial differential equations: A utility of fractional Wronskian. Math. Methods Appl. Sci. 2022, 45, 809–826. [Google Scholar] [CrossRef]
  9. He, J.H. Homotopy perturbation technique. Comput. Methods Appl. Mech. Eng. 1999, 178, 257–262. [Google Scholar] [CrossRef]
  10. He, J.H. Application of homotopy perturbation method to nonlinear wave equations. Chaos Solitons Fractals 2005, 26, 695–700. [Google Scholar] [CrossRef]
  11. He, J.H. Variational iteration method, a kind of nonlinear analytical technique. Some examples. Int. J. Non-Linear Mech. 1999, 34, 699–708. [Google Scholar] [CrossRef]
  12. Parsa, A.B.; Rashidi, M.M.; Beg, O.A.; Sardi, S.M. Semi-computational simulation of magneto-hemodynamic flow in a semi-porous channel using optimal homotopy and differential transform methods. Comput. Biol. Med. 2013, 43, 1142–1153. [Google Scholar] [CrossRef] [PubMed]
  13. Adomian, G.A. Review of the decomposition method in applied mathematics. J. Math. Anal. Appl. 1998, 135, 501–544. [Google Scholar] [CrossRef] [Green Version]
  14. Marinca, V.; Herisanu, N.; Bota, C. An optimal homotopy asymptotic method applied to the steady flow of a fourth-grade fluid past a porous plate. Appl. Math. Lett. 2009, 22, 245–251. [Google Scholar] [CrossRef] [Green Version]
  15. Srivastava, V.P. A Theoretical Model for Blood Flow in Small Vessels. Int. J. Appl. Appl. Math. 2007, 2, 51–65. [Google Scholar]
  16. Wang, C.Y.; Bassingthwaighte, J.B. Blood Flow in Small Curved Tubes. J. Biomech. Eng. 2003, 125, 910–913. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Bali, R.; Awasthi, U. Mathematical model of Blood Flow in Small Blood Vessel in the Presence of Magnetic Field. Appl. Math. 2001, 2, 264–269. [Google Scholar] [CrossRef] [Green Version]
  18. Yamamoto, T.; Nagayama, Y.; Tamura, M. A blood-oxygenation-dependent increase in blood viscosity due to a static magnetic field. Phys. Med. Biol. 2004, 49, 3267. [Google Scholar] [CrossRef] [PubMed]
  19. Alshare, A.; Tashtoush, B.; El-Khalil, H.H. Computational Modeling of Non-Newtonian Blood Flow Through Stenosed Arteries in the Presence of Magnetic Field. J. Biomech. Eng. 2013, 135, 1145–1153. [Google Scholar] [CrossRef]
  20. Weng, H.C. Hydrodynamic Modeling of Targeted Magnetic-Particle Delivery in a Blood Vessel. J. Biomech. Eng. 2013, 135, 034504. [Google Scholar] [CrossRef] [PubMed]
  21. Mekheimer, K.S. Peristaltic flow of blood under effect of a magnetic field in a non-uniform channels. Appl. Math. Comput. 2004, 153, 763–777. [Google Scholar] [CrossRef]
  22. Tenforde, T.S. Magnetically induced electric fields and currents in the circulatory system. Prog. Biophys. Mol. Biol. 2005, 87, 279–288. [Google Scholar] [CrossRef] [PubMed]
  23. Bhatti, M.M.; Zeeshan, A.; Bashir, F.; Sait, S.M.; Ellahi, R. Sinusoidal motion of small particles through a Darcy-Brinkman-Forchheimer microchannel filled with non-Newtonian fluid under electro-osmotic forces. J. Taibah Univ. Sci. 2021, 15, 514–529. [Google Scholar] [CrossRef]
  24. Rostami, S.; Ellahi, R.; Oztop, H.F.; Goldanlou, A.S. A study on the effect of magnetic field and the sinusoidal boundary condition on free convective heat transfer of non-Newtonian power-law fluid in a square enclosure with two constant-temperature obstacles using lattice Boltzmann method. J. Therm. Anal. Calorim. 2021, 144, 2557–2573. [Google Scholar] [CrossRef]
  25. Khan, A.S.; Xu, H.Y.; Khan, W. Magnetohydrodynamic Hybrid Nanofluid Flow Past an Exponentially Stretching Sheet with Slip Conditions. Mathematics 2021, 9, 3291. [Google Scholar] [CrossRef]
  26. Rehman, A.; Salleh, Z. Influence of Marangoni Convection on Magnetohydrodynamic Viscous Dissipation and Heat Transfer on Hybrid Nanofluids in a Rotating System among Two Surfaces. Mathematics 2021, 9, 2242. [Google Scholar] [CrossRef]
  27. Ali, B.; Naqvi, R.A.; Haider, A.; Hussain, D.; Hussain, S. Finite Element Study of MHD Impacts on the Rotating Flow of Casson Nanofluid with the Double Diffusion Cattaneo—Christov Heat Flux Model. Mathematics 2020, 8, 1555. [Google Scholar] [CrossRef]
  28. Ryu, J.; Hu, X.; Shadden, S.C. A Coupled Lumped-Parameter and Distributed Network Model for Cerebral Pulse-Wave Hemodynamics. J. Biomech. Eng. 2015, 137, 101009. [Google Scholar] [CrossRef] [Green Version]
  29. Srinivasacharya, D.; Rao, G.M. Mathematical model for blood flow through a bifurcated artery using couple stress fluid. Math. Biosci. 2016, 278, 37–47. [Google Scholar] [CrossRef]
  30. Sinha, A. MHD flow and heat transfer of a third order fluid in a porous channel with stretching wall: Application to hemodynamics. Alex. Eng. J. 2015, 54, 1243–1252. [Google Scholar] [CrossRef] [Green Version]
  31. Zaman, A.; Ali, N.; Sajid, M. Numerical simulation of pulsatile flow of blood in a porous-saturated overlapping stenosed artery. Math. Comput. Simul. 2017, 134, 1–16. [Google Scholar] [CrossRef]
  32. Caruntu, B.; Bota, C.; Bundau, O. Analytical simulation of magneto-hemodynamic flow in a semi-porous channel using the Polynomial Least Squares Method. ITM Web Conf. 2019, 29, 1–13. [Google Scholar] [CrossRef] [Green Version]
  33. Skalak, F.M.; Wang, C.Y. On the non-unique solutions of laminar flow through a porous tube or channel. SIAM J. Appl. Math. 1978, 34, 535–544. [Google Scholar] [CrossRef]
  34. Quaile, J.P.; Levy, E.K. Laminar flow in a porous tube with suction. Int. J. Heat Mass Transf. 1975, 97, 223–243. [Google Scholar] [CrossRef]
  35. Evans, E.A.; Skalak, R. Mechanics and Thermodynamics of Biomembranes. In Elsevier Biomedical; Hue, L., Van de Werve, G., Eds.; Elsevier Biomedical: Amsterdam, NY, USA, 1981; 464p. [Google Scholar]
Figure 1. Comparison of absolute errors corresponding to the approximation from [12] U D T M (red curve), [10] U H P M 3 terms (blue curve), [32] U P L S M (orange curve), and our LSHPM approximation U L S H P M (green curve).
Figure 1. Comparison of absolute errors corresponding to the approximation from [12] U D T M (red curve), [10] U H P M 3 terms (blue curve), [32] U P L S M (orange curve), and our LSHPM approximation U L S H P M (green curve).
Mathematics 10 00546 g001
Figure 2. Comparison of absolute errors corresponding to the approximation from [12] V D T M (red curve), [10] V H P M 3 terms (blue curve), [32] V P L S M (orange curve), and our LSHPM approximation V L S H P M (green curve).
Figure 2. Comparison of absolute errors corresponding to the approximation from [12] V D T M (red curve), [10] V H P M 3 terms (blue curve), [32] V P L S M (orange curve), and our LSHPM approximation V L S H P M (green curve).
Mathematics 10 00546 g002
Figure 3. The effect of the increase of R e on U for the case H a = 0 —three dimensional plot.
Figure 3. The effect of the increase of R e on U for the case H a = 0 —three dimensional plot.
Mathematics 10 00546 g003
Figure 4. The effect of the increase of R e on V for the case H a = 0 —three dimensional plot.
Figure 4. The effect of the increase of R e on V for the case H a = 0 —three dimensional plot.
Mathematics 10 00546 g004
Figure 5. The effect of the increase of H a on U for the case R e = 1 —three dimensional plot.
Figure 5. The effect of the increase of H a on U for the case R e = 1 —three dimensional plot.
Mathematics 10 00546 g005
Figure 6. The effect of the increase of H a on V for the case R e = 1 —three dimensional plot.
Figure 6. The effect of the increase of H a on V for the case R e = 1 —three dimensional plot.
Mathematics 10 00546 g006
Figure 7. The combined influence of R e and H a 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 7. The combined influence of R e and H a 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 .
Mathematics 10 00546 g007
Figure 8. The combined influence of R e and H a 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 .
Figure 8. The combined influence of R e and H a 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 .
Mathematics 10 00546 g008
Table 1. Comparison of the absolute errors of the approximate solutions U in case R e = 1 and H a = 0 .
Table 1. Comparison of the absolute errors of the approximate solutions U in case R e = 1 and H a = 0 .
y U HAM U DTM U PLSM U HPM 2 term . U HPM 3 t . U LSHPM 2 t . U LSHPM 3 t .
0.1 9.04 × 10 3 1.08 × 10 2 8.06 × 10 4 4.45 × 10 3 6.77 × 10 4 8.61 × 10 5 3.68 × 10 6
0.2 1.77 × 10 2 1.77 × 10 2 1.90 × 10 3 9.08 × 10 3 1.35 × 10 3 8.00 × 10 5 1.56 × 10 6
0.3 2.71 × 10 2 1.98 × 10 2 2.14 × 10 3 1.37 × 10 2 2.01 × 10 3 4.18 × 10 6 3.41 × 10 6
0.4 3.45 × 10 2 1.79 × 10 2 1.34 × 10 3 1.78 × 10 2 2.64 × 10 3 9.71 × 10 5 9.99 × 10 8
0.5 3.73 × 10 2 1.39 × 10 2 8.14 × 10 5 2.10 × 10 2 3.16 × 10 3 1.29 × 10 4 2.42 × 10 6
0.6 3.46 × 10 2 9.33 × 10 3 1.43 × 10 3 2.24 × 10 2 3.50 × 10 3 7.69 × 10 5 1.22 × 10 6
0.7 2.77 × 10 2 5.39 × 10 3 2.09 × 10 3 2.15 × 10 2 3.54 × 10 3 2.59 × 10 5 9.10 × 10 7
0.8 1.87 × 10 2 2.63 × 10 3 1.74 × 10 3 1.78 × 10 2 3.12 × 10 3 9.63 × 10 5 7.61 × 10 7
0.9 9.27 × 10 3 9.33 × 10 4 6.99 × 10 4 1.07 × 10 2 2.02 × 10 3 6.71 × 10 5 4.05 × 10 7
Table 2. Comparison of the absolute errors of the approximate solutions V in case R e = 1 and H a = 0 .
Table 2. Comparison of the absolute errors of the approximate solutions V in case R e = 1 and H a = 0 .
y V HAM V DTM V PLSM V HPM 2 term . V HPM 3 t . V LSHPM 2 t . V LSHPM 3 t .
0.1 2.21 × 10 4 1.19 × 10 6 1.48 × 10 7 2.25 × 10 5 2.13 × 10 6 1.41 × 10 6 2.20 × 10 8
0.2 5.75 × 10 4 4.10 × 10 6 7.17 × 10 7 1.14 × 10 4 6.55 × 10 6 1.21 × 10 5 9.78 × 10 9
0.3 6.47 × 10 2 7.91 × 10 6 1.34 × 10 6 2.83 × 10 4 9.89 × 10 6 3.20 × 10 5 4.69 × 10 9
0.4 3.03 × 10 4 1.17 × 10 5 1.26 × 10 6 4.95 × 10 4 9.25 × 10 6 5.07 × 10 5 1.40 × 10 8
0.5 2.63 × 10 4 1.49 × 10 5 4.45 × 10 7 6.84 × 10 4 3.56 × 10 6 5.66 × 10 5 1.53 × 10 8
0.6 7.02 × 10 4 1.65 × 10 5 3.80 × 10 7 7.75 × 10 4 5.54 × 10 6 4.58 × 10 5 1.84 × 10 8
0.7 7.62 × 10 4 1.59 × 10 5 5.63 × 10 7 7.12 × 10 4 1.37 × 10 5 2.53 × 10 5 1.67 × 10 8
0.8 4.73 × 10 4 1.23 × 10 5 2.24 × 10 7 4.87 × 10 4 1.54 × 10 5 7.48 × 10 6 1.46 × 10 8
0.9 1.25 × 10 4 5.88 × 10 6 2.11 × 10 8 1.80 × 10 4 8.05 × 10 6 2.77 × 10 7 1.41 × 10 8
Table 3. Comparison of the absolute errors of the approximate solutions for U for the case R e = 1 and H a = 1 .
Table 3. Comparison of the absolute errors of the approximate solutions for U for the case R e = 1 and H a = 1 .
y U HAM U DTM U PLSM U HPM 2 term . U HPM 3 t . U LSHPM 2 t . U LSHPM 3 t .
0.1 1.19 × 10 2 3.18 × 10 2 8.22 × 10 5 1.19 × 10 2 3.13 × 10 3 8.03 × 10 5 4.19 × 10 7
0.2 1.59 × 10 2 5.15 × 10 2 1.17 × 10 4 2.33 × 10 2 6.14 × 10 3 1.14 × 10 4 1.14 × 10 6
0.3 1.25 × 10 2 5.95 × 10 2 3.15 × 10 5 3.35 × 10 2 8.85 × 10 3 2.72 × 10 5 1.41 × 10 7
0.4 5.92 × 10 3 5.84 × 10 2 9.22 × 10 5 4.14 × 10 2 1.10 × 10 2 9.64 × 10 5 1.32 × 10 6
0.5 1.62 × 10 4 5.14 × 10 2 1.48 × 10 4 4.60 × 10 2 1.25 × 10 2 1.52 × 10 4 3.46 × 10 7
0.6 2.79 × 10 3 4.11 × 10 2 9.29 × 10 5 4.65 × 10 2 1.29 × 10 2 9.57 × 10 5 1.10 × 10 6
0.7 3.20 × 10 3 2.99 × 10 2 3.03 × 10 5 4.24 × 10 2 1.22 × 10 2 2.86 × 10 5 3.63 × 10 7
0.8 2.34 × 10 3 1.92 × 10 2 1.16 × 10 4 3.33 × 10 2 9.21 × 10 3 1.15 × 10 4 7.95 × 10 7
0.9 1.21 × 10 3 9.05 × 10 3 8.04 × 10 5 1.92 × 10 2 5.92 × 10 3 8.01 × 10 5 1.05 × 10 7
Table 4. Comparison of the absolute errors of the approximate solutions for V for the case R e = 1 and H a = 1 .
Table 4. Comparison of the absolute errors of the approximate solutions for V for the case R e = 1 and H a = 1 .
y V HAM V DTM V PLSM V HPM 2 term . V HPM 3 t . V LSHPM 2 t . V LSHPM 3 t .
0.1 4.90 × 10 4 2.70 × 10 4 3.06 × 10 8 7.74 × 10 5 1.53 × 10 6 1.14 × 10 8 1.91 × 10 8
0.2 1.15 × 10 3 5.83 × 10 4 8.27 × 10 8 3.03 × 10 4 8.01 × 10 6 5.24 × 10 6 7.27 × 10 9
0.3 1.19 × 10 3 5.51 × 10 4 1.37 × 10 7 6.22 × 10 4 2.06 × 10 6 1.87 × 10 5 4.18 × 10 10
0.4 4.36 × 10 4 1.45 × 10 4 1.04 × 10 7 9.38 × 10 4 3.79 × 10 6 3.38 × 10 5 3.16 × 10 9
0.5 7.25 × 10 4 4.47 × 10 4 1.51 × 10 8 1.14 × 10 3 5.51 × 10 6 4.05 × 10 5 1.35 × 10 9
0.6 1.66 × 10 3 9.52 × 10 4 9.11 × 10 8 1.16 × 10 3 6.56 × 10 6 3.41 × 10 5 3.91 × 10 10
0.7 1.90 × 10 3 1.12 × 10 3 6.89 × 10 8 9.70 × 10 4 6.30 × 10 5 1.91 × 10 5 1.24 × 10 9
0.8 1.35 × 10 3 8.72 × 10 4 1.55 × 10 8 6.05 × 10 4 4.48 × 10 5 5.55 × 10 6 7.28 × 10 10
0.9 4.73 × 10 4 3.42 × 10 4 6.71 × 10 10 2.04 × 10 4 1.71 × 10 6 1.06 × 10 7 3.35 × 10 9
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Paşca, M.S.; Bundău, O.; Juratoni, A.; Căruntu, 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

AMA Style

Paşca MS, Bundău O, Juratoni A, Căruntu B. The Least Squares Homotopy Perturbation Method for Systems of Differential Equations with Application to a Blood Flow Model. Mathematics. 2022; 10(4):546. https://doi.org/10.3390/math10040546

Chicago/Turabian Style

Paşca, Mădălina Sofia, Olivia Bundău, Adina Juratoni, and Bogdan Căruntu. 2022. "The Least Squares Homotopy Perturbation Method for Systems of Differential Equations with Application to a Blood Flow Model" Mathematics 10, no. 4: 546. https://doi.org/10.3390/math10040546

APA Style

Paşca, M. S., Bundău, O., Juratoni, A., & Căruntu, B. (2022). The Least Squares Homotopy Perturbation Method for Systems of Differential Equations with Application to a Blood Flow Model. Mathematics, 10(4), 546. https://doi.org/10.3390/math10040546

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop