Startseite Comparison of the method of variation of parameters to semi-analytical methods for solving nonlinear boundary value problems in engineering
Artikel Open Access

Comparison of the method of variation of parameters to semi-analytical methods for solving nonlinear boundary value problems in engineering

  • Travis J. Moore EMAIL logo und Vedat S. Ertürk
Veröffentlicht/Copyright: 19. Juli 2019
Veröffentlichen auch Sie bei De Gruyter Brill

Abstract

Solutions to nonlinear boundary value problems modelling physical phenomena in engineering applications have traditionally been approximated using numerical methods. More recently, several semi-analytical methods have been developed and used extensively in diverse engineering applications. This work compares the method of variation of parameters to semi-analytical methods for solving nonlinear boundary value problems arising in engineering. The accuracy and efficiency of the method of variation of parameters are compared to those of two widely used semi-analytical methods, the Adomian decomposition method and the differential transformation method, for three practical engineering boundary value problems: (1) the deflection of a cantilevered beam with a concentrated load, (2) an adiabatic tubular chemical reactor which processes an irreversible exothermal reaction, and (3) the electrohydrodynamic flow of a fluid in an ion drag configuration in a circular cylinder conduit. The accuracy and convergence of each method is investigated using the error remainder function. The method of variation of parameters is significantly more efficient than both the semi-analytical methods and traditional numerical methods while maintaining comparable accuracy. Unlike the semi-analytical methods, the efficiency of variation of parameters is independent of the nonlinearity. Variation of parameters is shown to be an attractive alternative to semi-analytical methods and traditional numerical methods for solving boundary value problems encountered in engineering applications in which solution efficiency is important.

1 Introduction

Mathematical modelling of physical phenomena in all branches of science and engineering frequently results in boundary value problems governed by nonlinear differential equations. These problems generally do not have closed-form, exact analytical solutions. However, because of their many important applications, significant efforts have been and continue to be made to accurately approximate solutions to these problems. Early methods developed to approximate solutions to nonlinear boundary value problems (BVPs) were strictly numerical. Whereas analytical methods are symbolic in nature, numerical techniques involve reformulating problems such that their solutions can be approximated with arithmetic operations whose error can be estimated [1]. Commonly used numerical methods for solving BVPs include shooting methods, Runge-Kutta methods, finite-difference methods, and integral equation methods. Extensive treatment of these and other numerical methods for solving BVPs are widely available [2, 3,].

More recently, several so-called “semi-analytical” methods for solving differential equations have been developed. Among these are the variational iteration method (VIM), the homotopy perturbation method (HPM), the homotopy analysis method (HAM), the Adomian decomposition method (ADM), and the differential transformation method (DTM). These methods have been used extensively in recent years and several texts and countless articles have been devoted to investigating their application to solve differential equations in numerous engineering applications from diverse fields such as solid mechanics, heat transfer, fluid mechanics, dynamics, and biomedicine [3, 4, 5, 6, 7, 8, 9, 10,]. The advantages of many of these semi-analytical approaches over numerical methods are their direct applicability to both linear and nonlinear equations without requiring linearization, discretization, or perturbation [11]. Additionally, they can be used to prove existence of solutions. However, these methods result in series solutions which must be truncated and can cause convergence problems. Although typically more accurate, semi-analytical methods are generally slower than numerical methods. This can be problematic in cases where solution efficiency is critical. An important example is inverse problems, in which the governing equation must be solved numerous times for optimization purposes [12].

A potential efficient and accurate alternative to semi-analytical and traditional numerical methods for solving nonlinear BVPs is the method of variation of parameters (VPM). This method was developed by Lagrange to solve linear, nonhomogeneous ordinary differential equations [13]. However, it can also be implemented to solve nonlinear differential equations in which the nonlinearity is found in the inhomogeneity [14, 15, 16, 17,]. The primary contribution of the present work is the comparison of the accuracy and efficiency of VPM to common semi-analytical methods. For reference, these methods will also be compared to traditional finite difference methods. It will be shown that VPM is a versatile and practical alternative for efficiently and accurately solving nonlinear BVPs that arise from models of physical phenomena in a wide variety of engineering applications.

Following an overview of how VPM, ADM, and DTM are implemented to solve nonlinear BVPs, solutions to four BVPs are presented using these methods in addition to traditional finite-difference methods. The first is a simple example with a known analytical solution so that the accuracy of each method can be quantified. Then, the solutions to three recently investigated engineering applications of nonlinear problems are shown. These include modelling (1) the deflection of a cantilevered beam with a concentrated load, (2) an adiabatic tubular chemical reactor which processes an irreversible exothermal reaction, and (3) the electrohydrodynamic flow of a fluid in an ion drag configuration in a circular cylinder conduit. These problems represent diverse fields of engineering whose solutions using semi-analytical methods have recently been published. Their governing equations have distinct linear and nonlinear portions which is used to demonstrate the versatility of the method of variation of parameters. The accuracy and efficiency of each method are compared for the three applications. Accuracy and convergence of each method is investigated using the error remainder function. The method of variation of parameters is shown to be an attractive alternative to semi-analytical methods and traditional numerical methods for solving boundary value problems encountered in engineering applications, especially for problems in which solution efficiency is an important consideration.

2 Method of variation of parameters

This section will demonstrate how the method of variation of parameters can be used to solve nonlinear differential equations in which the nonlinearity is in the inhomogeneity. For the sake of convenience and simplicity, the following description of the method is limited to second order differential equations. However, this method can be extended to higher order equations [18]. Consider the following equation

(1)d2ydx2+p(x)dydx+q(x)y=f(x,y,y),x0xx1

where f is a nonlinear function of the independent variable x aswell as the dependent variable y and/or its derivative y . The solution to this equation consists of the sum of the complementary and particular solutions. The complementary solution is the solution to the homogeneous equation corresponding to (1) that may be found using traditional ordinary, linear differential equation solution techniques [19]. By the principle of superposition, the complementary solution is yc = c1y1 + c2y2 where c1 and c2 are arbitrary constants and y1 and y2 form a fundamental set of solutions of the homogeneous equation. Herein lies the restriction that the homogeneous equation must be linear. dnydxn=fx,y,y,yn1 [18]. Therefore, this method can be used for any BVP of the form

The VPM method assumes that there exist functions u1 and u2 such that a particular solution has the form yp = u1y1 + u2y2 [19]. The unknown functions u1 and u2 are found by substituting the assumed particular solution into (1) while imposing the condition u1y1+u2y2=0,which results in a system of equations with unknowns u1 and u2.This system of equations can be solved to find expressions for u1and u2which can, in turn, be separated and integrated to yield the desired functions.

(2)u1=x0xft,y,yy2y1y2y2y1dt
(3)u2=x0xf(t,y,y)y1y1y2y2y1dt

The solution to (1) is therefore given by

(4)y(x)=c1y1+c2y2y1x0xf(t,y,y)y2y1y2y2y1dt+y2x0xf(t,y,y)y1y1y2y2y1dt

where the constants c1 and c2 are determined from the boundary conditions. Equation (4) is an exact solution to (1). However, due to the implicit nature of this solution, an iterative approach is required to approximate y. Numerical integration must be used during the iteration process and, if f is a function of y , finite difference methods are required to approximate the derivative. The accuracy of this method is therefore limited by the accuracy of the numerical methods employed. Details regarding the comparisons between exact solutions and those obtained by VPM are available in [20].

3 Adomian decomposition method

The Adomian decomposition method (ADM) is a common semi-analytical method used to solve nonlinear ordinary and partial differential equations [21]. It is applied by dividing the equation into linear and nonlinear parts and applying the inverse of the highest order linear derivative operator to both parts. The boundary conditions along with the source term are identified as the zeroth component of the series solution while the remaining terms are found by decomposing the nonlinear part into a series in which the terms are calculated recursively using Adomian polynomials. For the sake of consistency, the following description of this method is limited to second order equations of the form (1). However, the approach is applicable to nth order equations. Suppose that (1) is expressed as

(5)L(y)+R(y)=N(y)

where L=d2dx2is the second-order derivative operator, R = R=p(x)ddx+q(x)is the linear operator for the remainder of the linear portion of the equation, and N is the nonlinear operator N (y) = f (x, y, y). It is assumed that the inverse operator L−1 exists and can be taken as a definite integral of a function as follows.

(6)L1()=x0xx0t2()dt1dt2

Moving R (y) to the right side of (5) and applying the inverse operator to both sides yields

(7)y(x)=y0L1R(y)+L1N(y)

where y0 = y (x0) + y (x0) xy (x0) x0 is the first component of the decomposition series given by

(8)y(x)=n=0yn(x).

The nonlinear function N (y) is represented by the infinite series

(9)N(y)=n=0Nn

where the components Nn are the Adomian polynomials which are defined as [21]

(10)Nn(y0,y1,,yn)=1n!dndλn[N(i=0nyiλi)]λ=0,n0.

The remaining components of the series in (8) are determined recursively from the equation

(11)yn+1=L1Rn+L1Nn,n0.

ADM is a versatile method in that it can be used for linear and nonlinear ordinary and partial differential equations. However, it requires significant computational time to calculate the Adomian polynomials for higher order approximations. Additionally, this method is rapidly convergent in only very small regions while the convergence rate is much slower over wider regions [22].

4 Differential transformation method

The differential transformation method is another frequently used semi-analytical method. This method was first proposed by Zhou [23] to solve linear and nonlinear initial value problems arising in electric circuit analysis and has subsequently been used to solve numerous equations in a wide variety of applications [7, 24]. DTM is based on the Taylor series expansion and, like ADM, results in a series solution to the differential equation. Consider again the differential equation (1). An approximate solution to this equation is given by the series

(12)y(x)yN(x)=k=0NY(k)xk.

where the coefficients Y (0), Y (1)... Y (N) are obtained by applying the differential transform

(13)Y(k)=1k![dkydxk]x=x0

to both sides of the differential equation. The result is analgebraic equation whose unknowns are the coefficients of the series solution (12). Common operation properties of the one-dimensional differential transform can be found in [23, 24, 25,]. Equation (1) can be rearranged as follows

(14)d2ydx2=p(x)dydxq(x)y+f(x,y,y)

and the differential transform operations applied, resulting in the following recursive relation

(15)Y(k+2)=1(k+1)(k+2){l=0k[(l+1)Y(l+1)P(kl)+Y(l)Q(kl)]+F(k)}

where Y (k), P (k), Q (k), and F (k) are the differential transformations of y (x), p (x), q (x), and the nonlinear function f (x, y, y), respectively. This equation can be used along with the transformed boundary conditions to find the values of Y (k) required in the series solution. As is the case for ADM, DTM results in a truncated series solution which only provides an accurate approximation of the true solution in a very small region.

5 Example BVP: Comparison to analytical solution

The following example will demonstrate the use of all three methods to solve BVPs. Consider the following boundary value problem.

(16)d2ydx2+y=ex,y(0)=1,y(1)=0

The analytical solution is

(17)y(x)=12cos(x)+12(sin1ecos1)sin(x)+12ex.

If (16) is rearranged as

(18)d2ydx2=exy,y(0)=1,y(1)=0

it can be treated as an equation of the form of (1) with p (x) = q (x) = 0 and f = exy and can be solved using numerical and semi-analytical methods as presented below.

5.1 Solution by the method of variation of parameters

The complementary solution to (18) is yc = c1x + c2 where y1 = x and y2 = 1. Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c1 and c2 results in the following solution to (18).

(19)y(x)=1x01(ety)dt+x0x(ety)dt0xt(ety)dt

This solution is implicit and must be solved iteratively and the integrals must be evaluated numerically.

5.2 Solution by the Adomian decomposition method

Application of the inverse operator (6) to both sides of (18) yields

(20)y(x)=y0+xy0+L1(ex)L1(y)

where x0 = 0. Since y0 = y (0) = 1 and L−1 (ex) = exx−1, then from (11)

(21)yn+1=L1yn,n0.

The approximate solution to (18) by ADM is therefore given by

(22)y(x)=exx+y0xx22y0x36x5120+y0x5120+

where y0=y(0)is found by applying the boundary condition y (1) = 0.

5.3 Solution by the differential transformation method

Application of the differential transform operators to (18), where the differential transform of ex is 1/k!, gives the following recurrence relation.

(23)Y(k+2)=1(k+1)(k+2)(1k!Y(k))

From the differential transform given by (13), the transformed boundary conditions are

(24)Y(0)=1,k=0NkY(k)=0.

Utilizing the recurrence relation in (23) and the first transformed boundary condition, Y (k) for k ≥ 2 are found and the following series solution to (18) is obtained from (12).

(25)y(x)=1+Y(1)x+(16Y(1)6)x3+x424+Y(1)x5120+.

The second transformed boundary condition is used to solve for the unknown Y (1) = y (0).

5.4 Results and discussion

In order to properly compare the accuracy and efficiency of the method of variation of parameters to the semi-analytical methods, all three methods were implemented in MATLAB [26], a high-level scripting language. Equation (19) was solved iteratively with an initial guess of y (x) = 0 for 101 discrete values of x equally spaced over the domain. Simpson’s rules were used to perform the numerical integration. Convergence was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of y fell below 1 × 10−9. Further decreasing this criterion had no effect on the solution for all 15 significant figures of the double precision values. Convergence was achieved after 25 iterations which took 0.0338 seconds. The average percent error between the VPM solution and the exact solution at the discrete values was 3.89 × 10−7%.

To solve (18) in MATLAB using the Adomian decomposition method, the terms in the series solutions were solved for symbolically in terms of the unknown y0.At each iteration, the boundary condition y (1) = 0 was enforced to calculate a value for y0using a numeric solver and the series approximation of y was calculated over the domain at 101 discrete values of x. Convergence was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of yfell below 1 × 10−12. Further decreasing this criterion had no effect on the solution. Convergence was achieved in 3.24 seconds with 9 terms in the series solution. The average percent error between the ADM solution and the exact solution at the discrete values was 3.54 × 10−14%.

Implementation of the differential transformation method to solve (18) in MATLAB followed the same procedure as that used for ADM. The terms in the series were solved for in terms of the unknown constant Y (1), which was approximated at each iteration by applying the second boundary condition in (24). At each iteration, the series approximation of y was calculated over the domain at 101 discrete values of x. Convergence was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of y fell below 1 × 10−16. Further decreasing this criterion had no effect on the solution. Convergence was achieved in 0.047 seconds with 21 terms in the series solution. The average percent error between the DTM solution and the exact solution at the discrete values was 3.14 × 10−14%.

In order to compare these methods to traditional numerical methods, (18) was also solved with the “bvp4c” routine provided in MATLAB. This routine is a finite difference method (FDM) that implements the three-stage Lobatto IIIa collocation formula to solve boundary value problems. The solution to (18) was approximated at 101 discrete values of x. The relative tolerance was set to 1 × 10−7 below which there was no change to the solution. The computational time required was 0.278 seconds and the average percent error between the FDM solution and the exact solution at the discrete values was 1.35×10−7%. Table 1 shows a comparison of the computational time required to achieve convergence and the average percent error for each of the methods. It also compares the results of each method at six discrete points with 10 significant figures. Figure 1 shows a plot of the solutions at discrete points.

Fig. 1 Comparison of the solutions to equation (18) using all methods
Fig. 1

Comparison of the solutions to equation (18) using all methods

Table 1

Comparison of the required computational time to solve (18) and solutions at six data points

VPMADMDTMFDMExact
computational (seconds) time0.03383.240.04750.278NA
average error (%)3.89×10−73.54×10−143.14×10−141.35×10−7NA
y (0)1.00000000001.00000000001.00000000001.00000000001.0000000000
y (0.2)0.75568269400.75568269410.75568269410.75568269390.7556826941
y (0.4)0.53009503130.53009503130.53009503130.53009503090.5300950313
y (0.6)0.34304740550.34304740550.34304740550.34304740500.3430474055
y (0.8)0.21520863820.21520863820.21520863820.21520863760.2152086382
y (1)0.16781221440.16781221440.16781221440.16781221370.1678122144

All methods produced highly accurate results relative to the exact solution. The accuracies of the two semi-analytical methods were roughly the same and were much greater than those of the method of variation of parameters and the finite difference method, which also had comparable accuracies. VPM was the most efficient of all methods. The required computational time was much lower than ADM and FDM and slightly lower than DTM. The accuracy of VPM and FDM can be increased by using smaller discretization, which would decrease the efficiency. The accuracy of VPM can also be increased by using more accurate numerical quadrature methods. The following practical engineering applications show that the method of variation of parameters is, in general, significantly faster than semi-analytical methods and traditional numerical methods for soling boundary value problems while also maintaining a high level of accuracy.

6 Application 1: Deflection of a cantilevered beam with a concentrated load

Bisshopp and Drucker [27] obtained the following nonlinear boundary value problem as a model of the deflection of a cantilever beam under a concentrated load,

(26)d2θds2=λscosθ,θ(0)=0,θ(1)=0

where θ is the angle (radians) of the slope of the deflection curve at a point along the beam, s is the nondimensional arc length of the beam, and λ is the ratio of the load to the flexural rigidity of the beam (m−2). This model is applicable to beams that experience large deflections and is more realistic than those typically used in texts on mechanics of materials [28] because it does not neglect the square of the first derivative in the curvature formula when applying the Bernoulli-Euler theorem. Bisshopp and Drucker solved the equation using elliptical integrals [27].Na obtained numerical solutions using the finite difference method, the method of reduced physical parameters, and the method of invariant embedding [2]. More recently, a numerical solution was obtained by the collocation method with Haar wavelets [29]. The exact solution to (26) is unknown.

6.1 Solution by the method of variation of parameters

The complementary solution to (26) is θc = c1s + c2 where θ1 = s and θ2 = 1. Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c1 and c2 results in the following solution to (26),

(27)θ(s)=01f(t,θ)tdt01f(t,θ)dt+s0sf(t,θ)0sf(t,θ)tdt

where f (s, θ) = −λs cos θ is the nonlinear function.

6.2 Solution by the Adomian decomposition method

Application of the inverse operator to both sides of (26) and using (7) yields

(28)θ(s)=θ(0)+sθ(0)λsL1(cosθ)

where L−1 is given by (6) with s0 = 0. Since θ (0) = 0 then θ0 = θ (0) and from (11)

(29)θn+1=λsL1Nn,n0

where Nn (θ) = Nn (θ0, θ1, . . . , θn) are the Adomian polynomials of the nonlinear operator N (θ) = cos θ. The series solution can be found in terms of the unknown θ (0), which is determined by applying the boundary condition θ (1) = 0 to the approximate solution.

6.3 Solution by the differential transformation method

Taking the differential transform of both sides of (26) results in the following recursive relation

(30)Θ(k+2)=λ(k+1)(k+2)(l=0kδ(l1)G(kl))

where G (k) is the differential transform of cos θ which can be found in [30]. The boundary conditions are transformed as follows

(31)Θ(1)=0,k=0NΘ(k)=0.

The approximate series solution is obtained by substituting the expressions calculated from (30) into (12) and applying the first transformed boundary condition, resulting in an equation in terms of the unknown θ (0). This unknown is approximated from the second transformed boundary condition.

6.4 Results and discussion

Equation (26) was solved in MATLAB using each of the methods with λ = 8, consistent with [2, 27, 29]. The same procedures as those described in Section 5.4 were used to implement the four solution methods. Convergence for VPM was considered achieved when the Euclidean norm of the difference between the current and previous iteration values of θ fell below 1 × 10−13. Convergence was achieved after 103 iterations which took 0.0889 seconds. The convergence criterion for ADM was 1 × 10−5. The Adomian polynomials were calculated with the aid of the MATLAB function “allVl1” [31]. Convergence was achieved in 24.7 seconds with 11 terms in the series solution. Due to the extremely slow rate of convergence, the convergence criterion for DTM was set to 1×10−4. Convergence was achieved in 20.8 minutes with 28 terms in the series solution. Increasing the convergence criterion decreased the required computational time but also decreased the accuracy of the solution. The relative tolerance of the “bvp4c” routine for FDM was 1 × 10−9. The computational time required was 0.312 seconds. Table 2 shows a comparison of the computational time required to achieve convergence for each of the methods and compares the results at six discrete points with 10 significant figures. Figure 2 shows a plot of the solutions at discrete points.

Fig. 2 Comparison of the solutions to equation (26) using all methods
Fig. 2

Comparison of the solutions to equation (26) using all methods

Table 2

Comparison of the required computational time to solve (26) and solutions at six data points

VPMADMDTMFDM
computational time (seconds)0.088924.712480.312
θ (0)0.94012154010.94010280870.94012322420.9401215349
θ (0.2)0.93382068070.93380178720.93382238090.9338206748
θ (0.4)0.88910922690.88908919480.88911104410.8891092194
θ (0.6)0.76246172720.76243855850.76246386840.7624617169
θ (0.8)0.49452310890.49449407390.49452588710.4945230970
θ (1)0.00000000000.00000000000.00000000000.0000000000

Because the exact solution to (26) is unknown, the error remainder function is used as a measure of how well the approximate solutions satisfy the governing equation.

The error remainder function for (26) is

(32)ERn(s)=d2θnds2+λscosθn

where θn is the approximate solution to (26) after n iterations. The average error remainder function over the interval 0 ≤ s ≤ 1 at different values of n is a measure of the accuracy of the approximate solution as well as the rate of convergence. Figure 3 compares the average error remainder function as a function of iteration for the VPM, ADM, and DTM methods. Finite difference methods were employed to approximate θnin (32). The accuracy and rate of convergence for VPM is comparable to that of DTM while ADM has better accuracy and convergence.

Fig. 3 Comparison of the average error remainder as a function of iteration for the approximate solutions of equation (26) using VPM, ADM, and DTM.
Fig. 3

Comparison of the average error remainder as a function of iteration for the approximate solutions of equation (26) using VPM, ADM, and DTM.

7 Application 2: Adiabatic tubular chemical reactor

A mathematical model of the temperature distribution in an irreversible exothermic chemical reaction processed in an adiabatic tubular chemical reactor under steady state conditions is given by [32]

(33)d2ydx2λdydx=λμ(yβ)exp(y),y(0)=λy(0),y(1)=0

where y is the unknown nondimensional steady state temperature, x is the dimensionless length, λ is the Péclet number which represents the ratio of the rate of advection to the rate of diffusion, μ is the Damköhler number which represents the ratio of the chemical reaction rate to the rate of convective mass transport, and β is the dimensionless adiabatic temperature rise. The existence of the solution to this problem over various ranges has been investigated [33] and approximate solutions have been obtained using Adomian’s method for Hammerstein integral equations and the shooting method [34], the Sinc-Galerkin method [35], and the Sinc-collocation method [36].

7.1 Solution by the method of variation of parameters

The linear, homogeneous equation corresponding to (33) is solved by finding the roots of the auxiliary equation m2λm = 0 , resulting in the complementary solution yc = c1 + c2 exp (λx) where y1 = 1 and y2 = exp (λx). Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c1 and c2 results in the following solution to (33),

(34)y(x)=eλxλ[0xf(y)eλtdt01f(y)eλtdt]1λ0xf(y)dt

where f (y) = λμ (yβ) exp (y) is the nonlinear inhomogeneity.

7.2 Solution by the Adomian decomposition method

Application of the inverse operator to both sides of (33) and using (7) yields

(35)y(x)=y(0)+xy(0)+λ(L1y+μL1((yβ)exp(y)))

where L−1 is given by (6) with x0 = 0. The first boundary condition combined with the first two terms on the right side of (35) provide the first term in the series solution y0 = y (0) + λx y (0). The recursive relation is therefore

(36)yn+1=λ(L1yn+μL1(Nn)),n0

where Nn (y) = Nn (y0, y1, . . . , yn) are the Adomian polynomials of the nonlinear operator N (y) = (yβ) exp (y). The series solution can be found in terms of the unknown y (0), which is found by applying the second boundary condition to the approximate solution.

7.3 Solution by the differential transformation method

Taking the differential transform of both sides of (33) results in the following recursive relation

(37)Y(k+2)=1(k+1)(k+2)(λ(k+1)Y(k+1)+λμl=0k(Y(l)βδ(l)W(kl)))

where W (k) is the differential transform of exp (y) which can be found in [30]. The boundary conditions are transformed as follows

(38)Y(1)λY(0)=0,k=0NkY(k)=0.

The approximate series solution is obtained by substituting the expressions calculated from (37) into (12) and applying the first transformed boundary condition, resulting in an equation in terms of the unknown y (0). This unknown is approximated from the second transformed boundary condition.

7.4 Results and discussion

Equation (33) was solved in MATLAB using each of the methods with λ = 5, μ = 0.05, and β = 0.53, consistent with [34, 35, 36,]. The same procedures as those described in Section 5.4 were used to implement the four solution methods. The convergence criterion for VPM was 1×10−12 which was achieved in 0.0284 seconds after 8 iterations. Due to the slow rate of convergence of ADM, the convergence criterion set to 1 × 10−5. Convergence was achieved in 25.7 minutes with 6 terms in the series solution. The Adomian polynomials were calculated with the aid of the MATLAB function “allVl1” [31]. Increasing the convergence criterion decreased the required computational time but also decreased the accuracy of the solution. The convergence criterion for DTM was 1 × 10−7 which was achieved in 25.2 seconds with 20 terms in the series solution. An under-relaxation factor of 0.01 was required to achieve convergence. The relative tolerance of the “bvp4c” routine for FDM was 1 × 10−6. The computational time required was 0.225 seconds. Table 3 shows a comparison of the computational time required to achieve convergence for each of the methods and compares the results at six discrete points with 10 significant figures. Figure 4 shows a plot of the solutions at discrete points.

Fig. 4 Comparison of the solutions to equation (33) using all methods
Fig. 4

Comparison of the solutions to equation (33) using all methods

Table 3

Comparison of the required computational time to solve (33) and solutions at six data points

VPMADMDTMFDM
computational time (seconds)0.0284154125.20.225
y (0)0.0052161065200.0052020964030.0052161062160.005216106318
y (0.2)0.0103951499680.0103570015040.0103951494900.010395149767
y (0.4)0.0154481460430.0153446102400.0154481450870.015448145847
y (0.6)0.0202009424450.0199324965940.0202009401830.020200942263
y (0.8)0.0241782218850.0235797098820.0241782160690.024178221744
y (1)0.0260817693900.0251792722190.0260817553480.026081769381
Table 4

Comparison of the required computational time to solve (40) and solutions at six data points

VPMADMDTMFDM
computational time (seconds)0.02811.570.3980.381
w (0)0.0154428463590.0154433203600.0154433203020.015443317749
w (0.2)0.0148279127950.0148279129510.0148279128970.014828128286
w (0.4)0.0129805279790.0129805280510.0129805280090.012980622904
w (0.6)0.0098976775530.0098976775920.0098976775660.009897772768
w (0.8)0.0055735479780.0055735479950.0055735479830.005573765044
w (1)0.0000000000000.0000000000000.0000000000000.000000000000

The error remainder function for (33) is

(39)ERn(x)=d2yndx2λdyndxλμ(ynβ)exp(yn)

where yn is the approximate solution to (33) after n iterations. Figure 5 compares the average error remainder function over the interval 0 ≤ x ≤ 1 as a function of iteration for the VPM, ADM, and DTM methods. The accuracy and rate of convergence for ADM and DTM are comparable. The initial accuracy of VPM is better than the semi-analytical methods but converges after few iterations. This may be a result of the limitation of using finite-difference methods to approximate ynandynin (39) for VPM.

Fig. 5 Comparison of the average error remainder as a function of iteration for the approximate solutions of equation (33) using VPM, ADM, and DTM.
Fig. 5

Comparison of the average error remainder as a function of iteration for the approximate solutions of equation (33) using VPM, ADM, and DTM.

8 Application 3: Electrohydrodynamic flow

Electrohydrodynamics is the study of the dynamics of electrically charged fluids. Mathematical modelling of electrohydrodynamic flow is important in the design and analysis of flowmeters, accelerators, pumps, and magnetohydrodynamic generators. The electrohydrodynamic flow of a fluid in an “ion drag” configuration in a circular cylindrical conduit is governed by the nonlinear differential equation [37]

(40)d2wdr2+1rdwdr=H2(1w1αw),w(0)=0,w(1)=0

where w is the nondimensional fully-developed fluid velocity, r is the nondimensional radial distance from the center of the cylindrical conduit (the radius of the tube has been scaled to one), H is the Hartmann number which is the ratio of electromagnetic force to viscous force, and α is a dimensionless parameter related to the pressure gradient, the ion mobility, and the current density at the inlet of the conduit. Approximate solutions to (40) have been obtained using traditional numerical methods (Runge-Kutta/shooting method) [37] as well as with semi-analytical methods such as DTM [38] and the homotopy perturbation method [39]. The existence and uniqueness of the solution to (40) is examined in [40].

8.1 Solution by the method of variation of parameters

The linear, homogeneous equation corresponding to (40) is an Euler-Cauchy equation which can be solved by finding the roots of the auxiliary equation m2 = 0 , resulting in the complementary solution wc = c1 + c2 ln (r) where w1 = 1 and w2 = ln (r). Using the method of variation of parameters as outlined in Section 2 and applying the boundary conditions to find the constants c1 and c2 results in the following solution to (40),

(41)w(r)=01f(w)tln(t)dt0rf(w)tln(t)dt+ln(r)0rf(w)tdt

where f (w) = −H2 (1 − w/ (1 − αw)) is the nonlinear inhomogeneity.

8.2 Solution by the Adomian decomposition method

In (40) it is more efficient to keep the (1/r) dw/dr term on the left side of the equation such that the equation can be written as

(42)Lw=H2(1w1αw)

where the operator L is defined by

(43)L()=1rddr(rddr())

and the inverse operator is therefore

(44)L1()=1r1t20t2t1()dt1dt2.

Applying the inverse operator to both sides of (42) and enforcing the boundary conditions leads to the following recursive relation

(45)wn+1=L1(Nn),n0

where Nn (w) = Nn (w0, w1, . . . , wn) are the Adomian polynomials of the nonlinear operator N (w) = −H2 (1 − w/ (1 − αw)). By using the operator in (43) instead of the second-order derivative operator, application of the boundary conditions results in w0 = 0 thereby eliminating the need to use a solver to find w0. Equation (46) is used with (8) to determine the approximate series solution.

8.3 Solution by the differential transformation method

Taking the differential transform of both sides of (40) results in the following recursive relation

(46)W(k+1)=H2(k+1)2(δ(k1)F(k1))

where F (k) is the differential transform of 1 − w/ (1 − αw) which is given by

(47)F(k)={W(0)(1αW(0)),k=0w(k)l=0k1F(l)(δ(kl)αW(k1))1αW(0),k1.

The boundary conditions are transformed as follows

(48)W(1)=0,k=0NW(k)=0.

The approximate series solution to (40) is obtained by substituting the expressions calculated from (46) into (12) and applying the first transformed boundary condition, resulting in an equation in terms of the unknown w (0). This unknown is approximated from the second transformed boundary condition.

8.4 Results and discussion

Equation (40) was solved in MATLAB using each of the methods with H = 0.25 and α = 0.25, consistent with [37, 38, 39,]. The same procedures as those described in Section 5.4 were used to implement the four solution methods. The convergence criterion for VPM was 1 × 10−13 which was achieved in 0.0281 seconds after 8 iterations. The convergence criterion for ADM was 1 × 10−6 which was achieved in 1.57 seconds with 4 terms in the series solution. The Adomian polynomials were calculated with the aid of the MATLAB function “allVl1” [31]. The convergence criterion for DTM was 1 × 10−13 which was achieved in 0.398 seconds with 9 terms in the series solution. The relative tolerance of the “bvp4c” routine for FDM was 1 × 10−8. The computational time required was 0.381 seconds. Table 4 shows a comparison of the computational time required to achieve convergence for each of the methods and compares the results at six discrete points with 10 significant figures. Figure 6 shows a plot of the solutions at discrete points.

Fig. 6 Comparison of the solutions to equation (40) using all methods
Fig. 6

Comparison of the solutions to equation (40) using all methods

The error remainder function for (40) is

(49)ERn(r)=d2wndr2+1rdwndr=H2(1wn1αwn)

where wn is the approximate solution to (40) after n iterations. Figure 7 compares the average error remainder function over the interval 0 ≤ r ≤ 1 as a function of iteration for the VPM, ADM, and DTM methods. The accuracy and rate

Fig. 7 Comparison of the average error remainder as a function of iteration for the approximate solutions of equation (40) using VPM, ADM, and DTM.
Fig. 7

Comparison of the average error remainder as a function of iteration for the approximate solutions of equation (40) using VPM, ADM, and DTM.

of convergence for ADM and DTM are comparable. The initial accuracy of VPM comparable to that of the other methods but converges quickly while the accuracy of ADM and DTM continue to increase as the number of iterations increases. Again, this may be a result of the limitation of using finite-difference methods to approximate wnandwnin (49) for VPM.

9 Conclusion

This work has demonstrated the method of variation of parameters to be an accurate and very efficient method for solving nonlinear boundary value problems arising in engineering applications. The accuracy and efficiency of this method were compared to those of two semi-analytical methods, the Adomian decomposition method and the differential transformation method, as well as to traditional finite difference methods. Each of these methods were used to approximate solutions to BVPs encountered in diverse engineering applications whose governing equations had distinct linear and nonlinear portions. Generally, the accuracy of the semi-analytical methods was slightly better than that of the method of variation of parameters. At worst, the solutions from the different methods were equivalent up to at least five significant figures. However, the computational efficiency of the method of variation of parameters was greater than that of all other methods and was generally independent of the governing equation. In contrast, the efficiency of the semi-analytical methods was highly dependent on the nonlinearity of the BVP. The accuracy of the method of variation of parameters can be increased by using smaller discretization of the domain and by employing more accurate numerical quadrature methods to approximate the integrals. The method of variation of parameters is an attractive alternative to semi-analytical methods and traditional numerical methods for solving boundary value problems encountered in engineering applications, especially for problems in which solution efficiency is an important consideration.

References

[1] Chapra SC, Canale RP (2015) Numerical Methods for Engineers, 7th edn, McGraw Hill EducationSuche in Google Scholar

[2] Na TY (1979) Computational Methods in Engineering Boundary Value Problems. Academic PressSuche in Google Scholar

[3] Ardahaie SS, Amiri AJ, Amouei A, Hosseinzadeh Kh, Ganji DD (2018) Investigating the effect of adding nanoparticles to the blood flow in presence of magnetic field in a porous blood arterial. Informatics in Medicine Unlocked 10:71 – 8110.1016/j.imu.2017.10.007Suche in Google Scholar

[4] Ganji DD, Talarposhti RA (2017) Numerical and Analytical Solutions for Solving Nonlinear Equations in Heat Transfer. IGI Global10.4018/978-1-5225-2713-8Suche in Google Scholar

[5] Alizadeh M, Dogonchi AS, Gangi DD (2017) Nonlinear Heat Transfer in The Presence of Nanofluids. LAP Lambert Academic PublishingSuche in Google Scholar

[6] Hermann M, Saravi M (2016) Nonlinear Ordinary Differential Equations: Analytical Approximations and Numerical Methods. Springer10.1007/978-81-322-2812-7Suche in Google Scholar

[7] Hatami M, Ganji DD, Sheikholeslami M (2016) Differential Transformation Method for Mechanical Engineering Problems. Academic Press10.1016/B978-0-12-805190-0.00002-4Suche in Google Scholar

[8] Liao S (2014) Advances in the Homotopy Analysis Method. World Scientific Publishing Company10.1142/8939Suche in Google Scholar

[9] Haldar K (2015) Decomposition Analysis Method in Linear and Nonlinear Differential Equations. Chapman and Hall10.1201/b19396Suche in Google Scholar

[10] Barari A, Omidvar M, Ghotbi AR, Ganji DD (2008) Application of homotopy perturbation method and variational iteration method to nonlinear oscillator differential equations. Acta Applicandae Mathematicae 104(2):161 – 17110.1007/s10440-008-9248-9Suche in Google Scholar

[11] Moradi A, Ahmadikia H (2010) Analytical solution for different profiles of fin with temperature dependent thermal conductivity. Mathematical Problems in Engineering Article ID 568263:15 pages10.1155/2010/568263Suche in Google Scholar

[12] Müller JL, Siltanen S (2012) Linear and Nonlinear Inverse Problems with Practical Applications. SIAM10.1137/1.9781611972344Suche in Google Scholar

[13] Moulton FR (1984) An Introduction to Celestial Mechanics, 2nd revised edition. Dover PublicationsSuche in Google Scholar

[14] Bogdanova MP (1962) On a generalization of the method of variation of parameters. Doklady Akademii nauk BSSR 6:285 – 287Suche in Google Scholar

[15] Keckic J (1976) Additions to Kamkes treatise VII: Variation of parameters for nonlinear second order differential equations. Publikacije Elektrotehničkog fakulteta. Serija Matematika i fizika 544:31 – 36Suche in Google Scholar

[16] Mohyud-Din ST, Noor MA, Waheed A (2010) Variation of parameters method initial and boundary value problems. World Applied Sciences Journal 11(5):622 – 639Suche in Google Scholar

[17] Moore TJ, Jones MR (2015) Solving nonlinear heat transfer problems using variation of parameters. International Journal of Thermal Sciences 93:29 – 3510.1016/j.ijthermalsci.2015.02.002Suche in Google Scholar

[18] Noor MA, Mohyud-Din ST, Waheed A (2008) Variation of parameters method for solving fifth-order boundary value problems. Applied Mathematics and Information Sciences 2(2):135 – 141Suche in Google Scholar

[19] Zill DG, Cullen MR (2000) Advanced Engineering Mathematics, 2nd edn. Jones and Bartlett PublishersSuche in Google Scholar

[20] Moore TJ (2014) Application of variation of parameters to solve nonlinear multimode heat transfer problems. PhD dissertation, Brigham Young UniversitySuche in Google Scholar

[21] Adomian G (1994) Solving Frontier Problems of Physics: The Decomposition Method. Kluwer Academic10.1007/978-94-015-8289-6Suche in Google Scholar

[22] Urgulu Y, Kaya D, Inan IE (2011) Comparison of three semi-analytical methods for solving (1+1)-dimensional dispersive long wave equations. Computers and Mathematics with Applications 61:1278 – 129010.1016/j.camwa.2010.12.026Suche in Google Scholar

[23] Zhou JK (1986) Differential transformation and its application for electric circuits. Huazhong University Press, Wuhan, ChinaSuche in Google Scholar

[24] Bervillier C (2011) Status of the differential transformation method. Applied Mathematics and Computation 218:10158 – 1017010.1016/j.amc.2012.03.094Suche in Google Scholar

[25] Mirzaee F (2011) Differential transform method for solving linear and nonlinear systems of ordinary differential equations. Applied Mathematical Sciences 5:3465 – 3472Suche in Google Scholar

[26] MATLAB R2016b The MathWorks, Inc., Natick, Massachusetts, United StatesSuche in Google Scholar

[27] Bisshopp K, Drucker DC (1945) Large deflection of cantilever beams. Quarterly of Applied Mathematics 3(3):272 – 27510.1090/qam/13360Suche in Google Scholar

[28] Beer FP (2015) Mechanics of Materials, 7th edn. McGraw-Hill Higher EducationSuche in Google Scholar

[29] Siraj-ul-Islam, Aziz I, Šarler B (2010) The numerical solution of second-order boundary-value problems by collocation method with the Haar wavelets. Mathematical and Computer Modelling 52:1577 – 159010.1016/j.mcm.2010.06.023Suche in Google Scholar

[30] Chang S, Chang I (2008) A new algorithm for calculating one-dimensional differential transform of nonlinear functions. Applied Mathematics and Computation 195:799–80810.1016/j.amc.2007.05.026Suche in Google Scholar

[31] allVl1, Copyright © 2009, Bruno Luong, all rights reservedSuche in Google Scholar

[32] Poore AB (1989) A tubular chemical reactor model. A collection of nonlinear model problems contributed to the proceedings of the AMS-SIAM:28 – 31Suche in Google Scholar

[33] Feng W, Zhang G, Chai Y (2007) Existence of positive solutions for second order differential equations arising from chemical reactor theory. Discrete and Continuous Dynamical Systems:373 – 381Suche in Google Scholar

[34] Madbouly NM, McGhee DF, Roach GF (2001) Adomian’s method for Hammerstein integral equations arising from chemical reactor theory. Applied Mathematics and Computation 117:241 – 24910.1016/S0096-3003(99)00177-0Suche in Google Scholar

[35] Saadatmandi A, Razzaghi M, Dehghan M (2005) Sinc-Galerkin solution for nonlinear two point boundary value problems with applications to chemical reactor theory. Mathematical and Computer Modelling 42:1237 – 124410.1016/j.mcm.2005.04.008Suche in Google Scholar

[36] Rashidinia J, Nabati M (2013) Sinc-collocation solution for nonlinear two-point boundary value problems arising in chemical reactor theory. Malaya Journal of Matematik 4(1):97 – 106Suche in Google Scholar

[37] McKee S, Watson R, Cuminato JA, Caldwell J, Chen MS (1997) Calculation of electrohydrodynamic flow in a circular cylindrical conduit. Journal of Applied Mathematics and Mechanics 77:457 – 46510.1002/zamm.19970770612Suche in Google Scholar

[38] Mosayebidorcheh S (2013) Taylor series solution of the electrohydrodynamic flow equation. Journal of Mechanical Engineering and Technology 1(2):40 – 4510.18005/JMET0102001Suche in Google Scholar

[39] Khan NA, Jamil M, Mahmood A, Ara A (2012) Approximate Solution for the Electrohydrodynamic Flow in a Circular Cylindrical Conduit. ISRN Computational Mathematics, Volume 2012, Article ID 341069, 5 pages10.5402/2012/341069Suche in Google Scholar

[40] Paullet JE (1999) On the solutions of electrohydrodynamic flow in a circular cylindrical conduit. Journal of Applied Mathematics and Mechanics 79:357 – 36010.1002/(SICI)1521-4001(199905)79:5<357::AID-ZAMM357>3.0.CO;2-BSuche in Google Scholar

Received: 2018-10-07
Revised: 2019-03-10
Accepted: 2019-03-27
Published Online: 2019-07-19

© 2020 Travis J. Moore and Vedat S. Ertürk, published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.

Artikel in diesem Heft

  1. Comparison of the method of variation of parameters to semi-analytical methods for solving nonlinear boundary value problems in engineering
  2. Nonlinear H-infinity control for switched reluctance machines
  3. Energy flow of a 2018 FIA F1 racing car and proposed changes to the powertrain rules
  4. Risk index to monitor an anaerobic digester using a dynamic model based on dilution rate, temperature, and pH
  5. MHD Peristaltic flow of a nanofluid in a constricted artery for different shapes of nanosized particles
  6. Comparative study of homotopy perturbation transformation with homotopy perturbation Elzaki transform method for solving nonlinear fractional PDE
  7. Approximate method for solving strongly fractional nonlinear problems using fuzzy transform
  8. Numerical approach to MHD flow of power-law fluid on a stretching sheet with non-uniform heat source
  9. Entropy generation in an inclined porous channel with suction/injection
  10. Heat transfer from convecting-radiating fin through optimized Chebyshev polynomials with interior point algorithm
  11. Two dimensional simulation of laminar flow by three- jet in a semi-confined space
  12. Influence of temperature-dependent properties on a gravity-driven thin film along inclined plate
  13. Simulation and time-frequency analysis of the longitudinal train dynamics coupled with a nonlinear friction draft gear
  14. Study of differential transform technique for transient hydromagnetic Jeffrey fluid flow from a stretching sheet
  15. Generalized second-order slip for unsteady convective flow of a nanofluid: a utilization of Buongiorno’s two-component nonhomogeneous equilibrium model
  16. Numerical treatment for the solution of singularly perturbed pseudo-parabolic problem on an equidistributed grid
  17. Relative sea-level rise and land subsidence in Oceania from tide gauge and satellite GPS
  18. On finite series solutions of conformable time-fractional Cahn-Allen equation
  19. A generalized perspective of Fourier and Fick’s laws: Magnetized effects of Cattaneo-Christov models on transient nanofluid flow between two parallel plates with Brownian motion and thermophoresis
  20. MHD natural convection flow of Casson fluid in an annular microchannel containing porous medium with heat generation/absorption
  21. Numerical simulation of variable thermal conductivity on 3D flow of nanofluid over a stretching sheet
  22. Two meshless methods for solving nonlinear ordinary differential equations in engineering and applied sciences
  23. Thermoelastic analysis of FGM hollow cylinder for variable parameters and temperature distributions using FEM
  24. Qualitative analysis for two fractional difference equations
  25. MHD fractionalized Jeffrey fluid over an accelerated slipping porous plate
  26. Nonlinear analysis of high accuracy and reliability in traffic flow prediction
  27. Numerical solution of time-dependent Emden-Fowler equations using bivariate spectral collocation method on overlapping grids
  28. A reliable analytical technique for fractional Caudrey-Dodd-Gibbon equation with Mittag-Leffler kernel
  29. Accelerated HPSTM: An efficient semi-analytical technique for the solution of nonlinear PDE’s
  30. Effect of magnetized variable thermal conductivity on flow and heat transfer characteristics of unsteady Williamson fluid
  31. Couple stress fluid flow due to slow steady oscillations of a permeable sphere
  32. State-of-the-art of MW-level capacity oceanic current turbines
  33. Approximate solution for fractional attractor one-dimensional Keller-Segel equations using homotopy perturbation sumudu transform method
  34. Nonlinear absolute sea-level patterns in the long-term-trend tide gauges of the West Coast of North America
  35. Insight into the dynamics of non-Newtonian Casson fluid over a rotating non-uniform surface subject to Coriolis force
  36. Mixed convection flow in a vertical channel in the presence of wall conduction, variable thermal conductivity and viscosity
  37. A new structure formulations for cubic B-spline collocation method in three and four-dimensions
  38. Mathematical and numerical optimality of non-singular fractional approaches on free and forced linear oscillator
  39. MHD mixed convection on an inclined stretching plate in Darcy porous medium with Soret effect and variable surface conditions
  40. Comparative study of two techniques on some nonlinear problems based ussing conformable derivative
Heruntergeladen am 8.9.2025 von https://www.degruyterbrill.com/document/doi/10.1515/nleng-2018-0148/html
Button zum nach oben scrollen