Abstract
The Euler method is a typical one for numerically solving initial value problems of ordinary differential equations. Particle swarm optimization (PSO) is an efficient algorithm for obtaining the optimal solution of a nonlinear optimization problem. In this study, a PSO-based Euler-type method is proposed to solve the initial value problem of ordinary differential equations. In the typical Euler method, the equidistant grid points are always used to obtain the approximate solution. The existing shortcoming is that when the iteration number is increasing, the approximate solution could be greatly away from the exact one. Here, it is considered that the distribution of the grid nodes could affect the approximate solution of differential equations on the discrete points. The adopted grid points are assumed to be free and nonequidistant. An optimization problem is constructed and solved by particle swarm optimization (PSO) to determine the distribution of grid points. Through numerical computations, some comparisons are offered to reveal that the proposed method has great advantages and can overcome the existing shortcoming of the typical Euler formulae.
1. Introduction
Due to the complexity of differential equations arising in science and engineering, it is a continuously hot topic to propose various methods for numerically solving the initial and boundary value problems [1–7]. One of the techniques is to transform nonlinear differential equations to linear ordinary differential equations such as the homotopy analysis method [8]. The other of the techniques is to directly discrete differential equations [9]. It is worth noting that many typical methods have been proposed such as the finite difference and element ones [7, 10]. A basic technique of the two typical methods is to mesh the domain of the function governed by differential equations. Then, the approximate solutions on the nodes of the mesh are computed by a feasible approach. One can see that the mesh must be always predefined, which forms the base of the two-type numerical methods. In the present study, a typical initial value problem of the first order ordinary differential equation is considered and represented as follows:where the function is to be solved. is known, and is a constant. In order to numerically solve (1), it is worth noting that the most typical yet simple method is the Euler formula [11]. That is, under a series of nodes one has the following form:where the term denotes the approximate solution of In general, the nodes are always chosen to be equidistant with where h is the step size. Then, formula (2) can be further simplified as
It is noted that the Euler method has been extended and used widely [11–14]. One of the shortcomings is that the approximate solution obtained by the Euler method could be deviated greatly away from the exact one with a large iteration number unless the step size h is very small.
On the other hand, in many practical situations and the problems of sciences and engineering, various optimization problems have been created [15]. In order to obtain the optimal solutions, a great deal of typical methods have been proposed [15, 16] such as gradient methods, Newton’s method, conjugate direction methods, and others. Recently, the heuristic approaches have attracted much attention, such as genetic algorithms [17], particle swarm optimization (PSO) [18, 19], and others. In particular, it is worth noting that the PSO algorithm is efficient in solving many nonlinear optimization problems, and it has been developed widely [20–24] with application to various engineering problems [25–28]. The objective of the present paper is attempted to combine the typical Euler method and the PSO algorithm. Then, the existing shortcoming in the typical Euler method can be overcome. In addition, it is seen that the predefined mesh in the typical numerical methods could be not necessary. Hence, the mesh-free methods have been developed and applied to solve various engineering problems [4]. Inspired by the idea of mesh-free methods, we consider that the grid nodes are free and nonequidistant. When the Euler method is used, it is found that different distributions of the nodes lead to different values of There could exist an optimal distribution such that the approximate solution has the best accuracy. The remaining issue is how to find the optimal distribution of Motivated by the above considerations, here an optimization problem is constructed according to the error estimate. By using the particle swarm optimization [18], the constructed optimization problem is solved. Then, the nodes are obtained and used to determine the approximate solutions. Numerical results are carried out to show the advantages of the proposed approach through comparing with the typical Euler method. The observation reveals that a feasible distribution of the nodes can be indeed used to increase the accuracy of the approximate solution The paper is constructed as follows. Section 2 is focused on a brief review on the Euler method and its modified versions. In Section 3, for constructing an optimization problem, two approaches involving differential and integration techniques are provided. Section 4 offers numerical results by solving the optimization problems based on the PSO algorithm. By comparing with some known methods, the novelty and advantages of the proposed method are shown. The main conclusions are covered in Section 5.
2. A Brief Review on the Euler-Type Methods
Since the Euler formula (2) for initial value problems of ordinary differential equations is proposed, many modified ones have been addressed. It is convenient to call the Euler formula and its modified versions as the Euler-type methods. In what follows, let us offer some reviews on the Euler-type methods.
Now, we consider the case of with the equidistant nodes and give the following equality in terms of (1):
When in (4) is replaced by the Euler formula (3) can be obtained. Similarly, replacing by the implicit Euler formula can be determined as
It is seen that the implicit Euler formula (5) exhibits the same linear convergence as (3) with [11]. This means that formula (5) has no improvement with respect to (3) in accuracy. When applying the rectangular rule to the integration in (4), the trapezoidal rule is derived as
Moreover, the application of (3) to (6) leads to the improved Euler formula as follows:
It is found that formulae (6) and (7) have more accuracy than (3) and (5). In general, according to (4), the single-step method can be developed as follows:where is a real-valued function related to Hence, the single-step methods with higher order convergence can be constructed, such as the Runge–Kutta methods [11].
Furthermore, when we consider the case of with it follows
Then, the multistep methods can be developed as follows:where φ is a function given by Following the idea in (10), some explicit and implicit formulae have been constructed such as the Adams–Bashforth method, the Milne-Thomson method, and others [11]. The above Euler-type methods have been extended to solve fractional initial value problems [12], stochastic differential equations [13, 14], and others.
It is noted from the above-developed methods that they are based on the equidistant nodes to increase the accuracy of the typical Euler formula. Different from the above-mentioned methods, here we consider that the nodes are nonequidistant. By reasonably arranging the positions of the nodes, the accuracy of the typical Euler method can be increased. The proposed method can be further developed to improve the other numerical methods of differential and integral equations [7, 11].
3. The Methods of Constructing Optimization Problems
In this section, we construct several optimization problems to determine the grid points Two approaches are introduced in solving the initial value problem of ordinary differential equations (1).
3.1. Differential Method
First, the differential method is used to form an optimization problem. By considering the interval and the differential mean value theorem [29], there is at least one such that
It is seen that when the value of in (11) is approximated by using the Euler method (2) is given. When the value of in (11) is replaced by the implicit Euler method (5) is formed. On the other hand, if the value of can be determined explicitly, the application of (11) can lead to a perfect formula as follows:
Unfortunately, it is always difficult to give an explicit value of and which just forms the reason behind the Euler-type methods. On the other hand, we can change a way of thinking to construct an optimization problem as follows:
In other words, if the unknown is replaced by a constant the optimization problem (13) can be rewritten as
In what follows, let us analyze the optimization problem (14). It is convenient to assume meaning that the local truncation error is considered. Defining a new functionwe have the following result.
Theorem 1. Suppose that is a domain with If is continuously differentiated on the following estimate holds:where
Proof. According to (16), it followswhere and are located between and η lies between and Then, definingwe have the following relation:Under the consideration of the following result is determined:The proof is completed.
As shown in Theorem 1, it is found that the values of , and are dependent on the interval That is, a new function is constructed as follows:
In addition, considering the error estimate of the approximate solution we give the following result.
Theorem 2. Assume that the function is defined on with If is continuously differentiated on the error estimate is obtained as follows:
Proof. According to (14) and Theorem 1, we haveThe proof is completed.
Making use of Theorem 2, we can form a new optimization problem as follows:
Once the optimization problem (25) is solved, the nodes can be determined, and they are further used to determine through an Euler-type method such as one of (2)–(7). Hence, the approximate solution of is given aswhere
Of much interest is that the approximate solution could have more accuracy than the result obtained by the typical Euler method with the equidistant grid points. This means that when one wants to obtain an approximate solution with more accuracy, the optimization problem (25) can be solved to determine the distribution of grid points. The formats of the Euler-type methods (2)–(7) or their combinations can be used to give the approximate solution As compared with the typical Euler-type methods [11, 30], the main novelty of the present study is to assume that the grid points are free, and they are determined by constructing and solving an optimization problem.
3.2. Integration Method
Now, the integration method is introduced to construct an optimization problem. According to the initial value problem (1), we integrate both sides of the differential equation with respect to x from to The following Volterra integral equation is obtained:
By considering all the grid points integral equation (28) is rewritten as
Under the assumption of a continuous function on a domain there is at least one such thatwhere the integral mean value theorem has been applied [29]. When the values of and can be given explicitly, formula (30) is also explicit. However, it is always difficult to give the explicit values of and Then, it is assumed that and are replaced by and respectively. Hence, the optimization problem is constructed as follows:
Moreover, we introduce the Lipschitz conditions as follows:where and are positive Lipschitz constants. It gives
In addition, applying the differential mean value theorem, there is at least one such thatwhere is located between x and By consideringwe have
Applying the Cauchy–Schwarz inequality to (34) together with (33), the following relation is given:
Then, the optimization problem (31) is transformed to the following form:
It is easy to give the optimal solution of (38) aswith This means that when the optimization problem (38) is solved, the approximate solution of can be determined as follows:
In other words, midpoint formula (40) has been retrieved by solving an optimization problem.
Furthermore, of much important is that here the distribution of the grid points is considered to be free. Hence, a new optimization problem is constructed as follows:
Once optimization problem (41) is solved, the distribution of the grid point is determined. Formulae (2)–(7) and (40) can be used to obtain the approximation solution . It is hopeful to give an approximate solution with more accuracy than that by using the equidistant grid points. Moreover, it is seen from (40) that there is an unknown term For the sake of simplicity, the Euler formula is applied to give
Formula (40) is rewritten as
As compared with the differential method, the obtained functions (25) and (39) are similar, when the Lipschitz constants and are replaced by and , respectively.
According to the above analysis, the two optimization problems (25) and (41) are constructed to determine the nonequivalent grids such that the approximate solution to has more accuracy than that by using the equivalent grids. Furthermore, when we consider the grids and the error estimate on all the grids the optimization problem (13) can be reread as
In terms of (44) and (25), the derived optimization problem is obtained as
In the following computations, the optimization problems (25) and (45) are all used to determine the grid points. Some comparisons are offered to illustrate the advantages of the proposed methods in giving a good approximation solution.
4. Numerical Results and Discussion
In order to show the effectiveness of the proposed method, some numerical results are reported in this section.
4.1. The Initial Value Problems of the First Order Ordinary Differential Equations
Example 1. Consider the following initial value problem of the first order differential equation [30]:with the exact solutionFirst, let us construct an optimization problem according to (25). In terms of (46), it givesmeaning thatWe further have the following estimates:Then, it followsOne can find from (51) that the value of is unknown. Hence, the main term of (51) related to the nodes is extracted to yield the following function:The optimization problemis used to determine the values of the nodes. Hence, the application of (43) leads to the computation formula as follows:Second, let us focus on the optimal solution of (53). The particle swarm optimization developed by Kennedy and Eberhart [18, 19] is utilized to achieve the objective. For convenience, the formulae of updating the position and velocity of a particle are given as follows [31]:The meanings of the terms in the above formulae are explained in Table 1.
Moreover, for numerical computations, it is assumed that ω is monotonically decreasing from 0.9 to 0.4 with respect to the generation number and When applying the PSO algorithm to the optimization problem (53), a particle is considered as the following -dimensional vector:subject to the constraint condition:where c is a prescribed constant. The size of particle swarm is 1000, and the maximum of the generation is 200. The symbols and stand for the individual best position of a particle and the global best position of 1000 particles, respectively.
Figure 1 is drawn to show the variation of the values of F versus the generation number with and It is seen from Figure 1 that the values of F decrease rapidly and then tend to a constant with the increasing generation number. This means that a stable and optimal solution has been obtained when the generation number is 200. The computed nonequidistant grid points are given as follows:The approximate solution using the nonequidistant grid points and formula (54) is obtained as When the equidistant grid points and (2) are used, we obtain The exact solution is Clearly, the solution is more approximate to than The observation reveals that the proposed Euler-type method with optimization is effective to obtain an approximate solution with more accuracy. In addition, under the condition we also choose for computations, and the observations are shown in Table 2. It is found from Table 2 that the approximate solutions of determined by the Euler method are away from the exact one when the values of are increasing. This is also attributed to the fact that the number of the chosen grid points is only 6. Fortunately, the proposed method is still effective regardless of the increasing values of In other words, it is interesting to find that the obtained results by using the proposed method are in good approximation to the exact ones. This means that the shortcoming of the typical Euler method can be overcome by using the developed method in the present study.
In what follows, we construct the other optimization problem in terms of (45) as follows:withFor the sake of comparisons, we choose and for numerical computations. After running the PSO algorithm, the optimal solution is obtained, and the nonequivalent grids are given asThe obtained approximation solution for is Expectedly, it has more errors than the value in Table 2 by solving the optimization problem (53). For convenience, Figure 2 is drawn to show the variations of the absolute errors and versus the discrete grids, where and stand for the approximate solutions by solving (53) and (54), respectively. It is seen from Figure 2 that the absolute errors are increasing with the increasing values of the discrete points. The optimization problem (53) focuses on the minimization of the error on the point The aim of the optimization problem (54) is the minimization of the global error on all discrete points.
Example 2. In engineering applications, the widely studied Riccati differential equation is expressed as the following form [32, 33]:Here, we assume that it is subject to the initial value condition with , and Then, the explicit solution is given as [34, 35]Under the consideration of (46), we haveIn virtue of (25), an optimization problem is constructed as follows:For numerical computations, we choose with By running the PSO algorithm, the optimization problem (65) can be solved, and the grid points are determined. Then, by using (43), the approximate solutions are obtained asIn order to show the advantages of the proposed approach, the approximate solutions in [34, 35] are recalled and given in the following:As shown in Table 3, the absolute errors , and are computed by virtue of different approaches. It is found from Table 3 that the better approximation solutions have been given by using the proposed method as compared to the known ones [34, 35]. In particular, the obtained results are stable, and the absolute errors are controlled in an acceptable domain.
4.2. The Initial Value Problem of Higher Order Ordinary Differential Equations
It is obvious that the developed method in the previous section can be used to numerically solve the initial value problem of the higher order differential equations. For example, as shown in [36], we consider a damped oscillator circuit for reflecting the behaviour of a stiffness system.
Example 3. The initial value problem of the second order ordinary differential equation is expressed as follows:where I is the electric current; and C are constants with and respectively. The explicit solution can be computed asIn order to apply the proposed method, it is assumed thatThe initial problem (68) is rewritten as the system of the first order differential equations: subject to the initial value conditions:In addition, we can also reexpress (71) as the style of the vector. That is, it is supposed thatThen, one haswhereFor a series of grid points the application of (43) leads towhere denotes the approximate solution of When the optimal solution is focused on, the grid points are determined by solving the following optimization problem:When considering the global optimal solution, the optimization problem for determining the grid points is constructed as follows:In what follows, the grid points are chosen for numerical computations. By running the PSO algorithm, the optimal solutions of (77) and (78) are obtained, and the following grid points are determined asrespectively. One can see that the optimal solution of (77) yields the equivalent grid points. Figure 3 shows the variations of the approximate and explicit solutions versus the time The obtained results reveal that the approximate solutions are away from the explicit one with the increasing time. When the proposed methods are applied, the approximate solutions have more accuracy than those by using the typical Euler method.
5. Conclusions
In the typical Euler-type methods for numerically solving ordinary differential equations, the grid points are always assumed to be equidistant. When the iteration number is increasing, the approximate solution could be greatly away from the exact one. Here, we have improved the typical Euler method by reasonably arranging the nonequidistant nodes. Numerical results have been carried out to verify the novelty and advantages of the developed method by comparing with the known ones. The main findings are given as follows:(i)A feasible distribution of the nonequidistant nodes can increase the accuracy of the typical Euler method based on the equidistant ones(ii)The positions of the nonequidistant nodes can be determined by constructing and solving an optimization problem related to the error estimate(iii)The shortcoming of the typical Euler method can be overcome to a certain degree
It is seen that the developed method is based on the novel idea different from those in the known Euler-type methods. In the future works, the idea in the proposed method could be used to develop the finite element method and the difference method in numerically solving ordinary and partial differential equations. The known software for numerically solving engineering and scientific problems could be further improved with more accuracy and efficiency.
Data Availability
The data used to support the findings of this study are included in this paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work was supported by the National Natural Science Foundation of China (nos. 11872155 and 11362002), the Guangxi Natural Science Foundation (no. 2016GXNSFAA380261), the 2017 Guangxi High School Innovation Team and Outstanding Scholars Plan, and the Innovation Project of Guangxi Graduate Education (no. YCSW2019045).