# A Particle Swarm Optimization-Based Method for Numerically Solving Ordinary Differential Equations

Authors: Xian-Ci Zhong; Jia-Ye Chen; Zhou-Yang Fan
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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:(1)duxdx=fx,u,x≥x0,ux0=u0,where the function ux is to be solved. fx,u is known, and u0 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 x0<x1<x2<⋯<xk<⋯, one has the following form:(2)uk+1=uk+xk+1−xkfxk,uk,k=0,1,2,…,where the term uk denotes the approximate solution of uxk. In general, the nodes are always chosen to be equidistant with xk=x0+kh, where h is the step size. Then, formula (2) can be further simplified as(3)uk+1=uk+hfxk,uk,k=0,1,2,….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 k, 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 x0,x1,…,xk lead to different values of uk+1. There could exist an optimal distribution such that the approximate solution uk+1 has the best accuracy. The remaining issue is how to find the optimal distribution of x0,x1,…,xk. 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 x0,x1,…,xk can be indeed used to increase the accuracy of the approximate solution uk+1. 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 ofx∈xk,xk+1 with the equidistant nodes xk=x0+kh and give the following equality in terms of (1):(4)uxk+1=uxk+∫xkxk+1fx,udx.Whenfx,u in (4) is replaced by fxk,uk, the Euler formula (3) can be obtained. Similarly, replacing fx,u by fxk+1,uk+1, the implicit Euler formula can be determined as(5)uk+1=uk+hfxk+1,uk+1,k=0,1,2,….It is seen that the implicit Euler formula (5) exhibits the same linear convergence as (3) with h⟶0 [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(6)uk+1=uk+h2fxk,uk+fxk+1,uk+1,k=0,1,2,….Moreover, the application of (3) to (6) leads to the improved Euler formula as follows:(7)uk+1=uk+h2fxk,uk+fxk+1,uk+hfxk,uk,k=0,1,2,….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:(8)uk+1=uk+hφxk,uk;h,k=0,1,2,…,where φxk,uk;h is a real-valued function related to fx,u. 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 ofx∈xk+r−j,xk+r with 1≤j≤r, it follows(9)uxk+r=uxk+r−j+∫xk+r−jxk+rfx,udx.Then, the multistep methods can be developed as follows:(10)uk+r+ar−1uk+r−1+⋯+a0uk=hφxk,uk,…,uk+r−1;h,where φ is a function given by f. 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 pointsx0<x1<x2<⋯<xn. 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 intervalxk,xk+1 and the differential mean value theorem [29], there is at least one ξk∈xk,xk+1 such that(11)uxk+1−uxkxk+1−xk=dudxx=ξk=fξk,uξk,k=0,1,2,….It is seen that when the value ofξk in (11) is approximated by using xk, the Euler method (2) is given. When the value of ξk in (11) is replaced by xk+1, the implicit Euler method (5) is formed. On the other hand, if the value of ξk can be determined explicitly, the application of (11) can lead to a perfect formula as follows:(12)uk+1=uk+xk+1−xkfξk,uξk.Unfortunately, it is always difficult to give an explicit value ofξk and uξk, 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:(13)minuxk+1−uk+1.In other words, if the unknownξk is replaced by a constant ck∈xk,xk+1, the optimization problem (13) can be rewritten as(14)minuxk−uk+xk+1−xkfξk,uξk−fck,uck.In what follows, let us analyze the optimization problem (14). It is convenient to assume uxk=uk, meaning that the local truncation error is considered. Defining a new function(15)Qxk,xk+1,ξk,ck=xk+1−xkfξk,uξk−fck,uck,we have the following result.Theorem 1.
Suppose thatD⊂R2 is a domain with fx,u:D⟶R. If fx,u is continuously differentiated on D, the following estimate holds:(16)Qxk,xk+1,ξk,ck≤xk+1−xk2Mk+Nk⋅Lk,where(17)Mk=maxx∈xk,xk+1∂f∂x,Lk=maxx∈xk,xk+1f,Nk=maxx∈xk,xk+1,x,u∈D∂f∂u.Proof.
According to (16), it follows(18)fξk,uξk−fck,uck=fξk,uξk−fck,uξk+fck,uξk−fck,uck=fξk,uξk−fck,uξk+fck,uξk−fck,uck=∂fξ¯k,uξk∂xξk−ck+∂fck,η∂uuξk−uck=∂fξ¯k,uξk∂xξk−ck+∂fck,η∂u⋅fηk,uηk⋅ξk−ck,where ξ¯k and ηk are located between ξk and ck; η lies between uξk and uck. Then, defining(19)Mk=maxx∈xk,xk+1∂f∂x,Lk=maxx∈xk,xk+1f,Nk=maxx∈xk,xk+1,x,u∈D∂f∂u,we have the following relation:(20)fξk,uξk−fck,uck≤ξk−ckMk+Nk⋅Lk.
Under the consideration ofξk,ck∈xk,xk+1, the following result is determined:(21)Qxk,xk+1,ξk,ck≤xk+1−xk2Mk+Nk⋅Lk.
The proof is completed.As shown in Theorem 1, it is found that the values ofMk, Nk, and Lk are dependent on the interval xk,xk+1. That is, a new function is constructed as follows:(22)Q¯kxk,xk+1=xk+1−xk2Mk+Nk⋅Lk.In addition, considering the error estimate of the approximate solutionuk+1, we give the following result.Theorem 2.
Assume that the functionfx,u is defined on D⊂R2 with fx,u:D⟶R. If fx,u is continuously differentiated on D, the error estimate is obtained as follows:(23)uxk+1−uk+1≤∑j=0kQ¯jxj,xj+1.Proof.
According to (14) and Theorem 1, we have(24)uxk+1−uk+1=uxk−uk+xk+1−xkfξk,uξk−fck,uck≤uxk−uk+Q¯kxk,xk+1≤uxk−1−uk−1+Q¯kxk,xk+1+Q¯k−1xk−1,xk⋯⋯⋯⋯⋯≤ux0−u0+Q¯kxk,xk+1+Q¯k−1xk−1,xk+⋯+Q¯0x0,x1=∑j=0kQ¯jxj,xj+1.
The proof is completed.Making use of Theorem 2, we can form a new optimization problem as follows:(25)minx0<x1<⋯<xk+1∑j=0kQ¯jxj,xj+1.Once the optimization problem (25) is solved, the nodes xjj=1,2,…,k+1 can be determined, and they are further used to determine uj through an Euler-type method such as one of (2)–(7). Hence, the approximate solution of ux x0≤x≤xk+1 is given as(26)uax=∑j=0k+1ujljx,where(27)l0x=x−x1x0−x1,x∈x0,x1,0,x∈x1,xk+1,ljx=x−xj−1xj−xj−1,x∈xj−1,xj,x−xj+1xj−xj+1,x∈xj,xj+1,0,x∉xj−1,xj+1,lk+1x=x−xkxk+1−xk,x∈xk,xk+1,0,x∈x0,xk.Of much interest is that the approximate solutionuk+1 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 uk+1 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 uk+1. 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 x0 to xk+1. The following Volterra integral equation is obtained:(28)uxk+1=u0+∫x0xk+1fx,udx.By considering all the grid pointsx0<x1<x2<⋯<xk+1, integral equation (28) is rewritten as(29)uxk+1=u0+∑j=0k∫xjxj+1fx,udx.Under the assumption of a continuous functionfx,u on a domain D, there is at least one ξj∈xj,xj+1 such that(30)uxk+1=u0+∑j=0kxj+1−xjfξj,uξj,where the integral mean value theorem has been applied [29]. When the values of ξj and uξj can be given explicitly, formula (30) is also explicit. However, it is always difficult to give the explicit values of ξj and uξj. Then, it is assumed that ξj and uξj are replaced by cj and ucj, respectively. Hence, the optimization problem is constructed as follows:(31)minuxk+1−uk+1=mincj∈xj,xj+1∑j=0kxj+1−xjfξj,uξj−fcj,ucj=mincj∈xj,xj+1∑j=0k∫xjxj+1fx,udx−∑j=0kxj+1−xjfcj,ucj=mincj∈xj,xj+1∑j=0k∫xjxj+1fx,udx−∑j=0k∫xjxj+1fcj,ucjdx=mincj∈xj,xj+1∑j=0k∫xjxj+1fx,u−fcj,ucjdx.Moreover, we introduce the Lipschitz conditions as follows:(32)fx,u−fcj,u≤L¯jx−cj,x,cj∈xj,xj+1,∀u,fx,ux−fx,ucj≤L˜jux−ucj,∀x∈xj,xj+1,where L¯j and L˜j are positive Lipschitz constants. It gives (33)∑j=0k∫xjxj+1fx,u−fcj,ucjdx=∑j=0k∫xjxj+1fx,u−fcj,ux+fcj,ux−fcj,ucjdx≤∑j=0k∫xjxj+1fx,u−fcj,ux+fcj,ux−fcj,ucjdx=∑j=0k∫xjxj+1L¯jx−cj+L˜jux−ucjdx.In addition, applying the differential mean value theorem, there is at least oneηj such that(34)ux−ucj=fηj,uηj⋅x−cj,where ηj is located between x and cj. By considering(35)Lj=maxx∈xj,xj+1f,we have(36)∑j=0k∫xjxj+1L¯jx−cj+L˜jux−ucjdx=∑j=0k∫xjxj+1L¯jx−cj+L˜jfηj,uηj⋅x−cjdx≤∑j=0k∫xjxj+1L¯j+Lj⋅L˜jx−cjdx.Applying the Cauchy–Schwarz inequality to (34) together with (33), the following relation is given:(37)∑j=0k∫xjxj+1fx,u−fcj,ucjdx≤∑j=0kL¯j+Lj⋅L˜j∫xjxj+1x−cjdx≤∑j=0kL¯j+Lj⋅L˜j∫xjxj+112dx⋅∫xjxj+1x−cj2dx1/2=∑j=0kL¯j+Lj⋅L˜jxj+1−xj3xj+1−cj3+xj−cj31/2.Then, the optimization problem (31) is transformed to the following form:(38)mincj∈xj,xj+1∑j=0kL¯j+Lj⋅L˜jxj+1−xj3xj+1−cj3+xj−cj31/2.It is easy to give the optimal solution of (38) as(39)Φ=36∑j=0kL¯j+Lj⋅L˜jxj+1−xj2,with cj=xj+xj+1/2. This means that when the optimization problem (38) is solved, the approximate solution of uxk+1 can be determined as follows:(40)uk+1=u0+∑j=0kxj+1−xjfxj+xj+12,uxj+xj+12.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:(41)F=minx0<x1<⋯<xk+1Φ.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 uk+1. 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 uxj+xj+1/2. For the sake of simplicity, the Euler formula is applied to give(42)uxj+xj+12≈uk+xj+1−xj2fxk,uk.Formula (40) is rewritten as(43)uk+1=u0+∑j=0kxj+1−xjfxj+xj+12,uj+xj+1−xj2fxj,uj.As compared with the differential method, the obtained functions (25) and (39) are similar, when the Lipschitz constants L¯j and L˜j are replaced by Mk and Nk, respectively.According to the above analysis, the two optimization problems (25) and (41) are constructed to determine the nonequivalent grids x0<x1<⋯<xk such that the approximate solution uk+1 to uxk+1 has more accuracy than that by using the equivalent grids. Furthermore, when we consider the grids a=x0<x1<⋯<xn=b and the error estimate on all the grids xk k=1,2,…,n, the optimization problem (13) can be reread as(44)min∑k=0n−1uxk+1−uk+1.In terms of (44) and (25), the derived optimization problem is obtained as(45)minx0<x1<⋯<xn∑k=0n−1∑j=0kQ¯jxj,xj+1.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.


## 3.1. Differential Method

First, the differential method is used to form an optimization problem. By considering the intervalxk,xk+1 and the differential mean value theorem [29], there is at least one ξk∈xk,xk+1 such that(11)uxk+1−uxkxk+1−xk=dudxx=ξk=fξk,uξk,k=0,1,2,….It is seen that when the value ofξk in (11) is approximated by using xk, the Euler method (2) is given. When the value of ξk in (11) is replaced by xk+1, the implicit Euler method (5) is formed. On the other hand, if the value of ξk can be determined explicitly, the application of (11) can lead to a perfect formula as follows:(12)uk+1=uk+xk+1−xkfξk,uξk.Unfortunately, it is always difficult to give an explicit value ofξk and uξk, 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:(13)minuxk+1−uk+1.In other words, if the unknownξk is replaced by a constant ck∈xk,xk+1, the optimization problem (13) can be rewritten as(14)minuxk−uk+xk+1−xkfξk,uξk−fck,uck.In what follows, let us analyze the optimization problem (14). It is convenient to assume uxk=uk, meaning that the local truncation error is considered. Defining a new function(15)Qxk,xk+1,ξk,ck=xk+1−xkfξk,uξk−fck,uck,we have the following result.Theorem 1.
Suppose thatD⊂R2 is a domain with fx,u:D⟶R. If fx,u is continuously differentiated on D, the following estimate holds:(16)Qxk,xk+1,ξk,ck≤xk+1−xk2Mk+Nk⋅Lk,where(17)Mk=maxx∈xk,xk+1∂f∂x,Lk=maxx∈xk,xk+1f,Nk=maxx∈xk,xk+1,x,u∈D∂f∂u.Proof.
According to (16), it follows(18)fξk,uξk−fck,uck=fξk,uξk−fck,uξk+fck,uξk−fck,uck=fξk,uξk−fck,uξk+fck,uξk−fck,uck=∂fξ¯k,uξk∂xξk−ck+∂fck,η∂uuξk−uck=∂fξ¯k,uξk∂xξk−ck+∂fck,η∂u⋅fηk,uηk⋅ξk−ck,where ξ¯k and ηk are located between ξk and ck; η lies between uξk and uck. Then, defining(19)Mk=maxx∈xk,xk+1∂f∂x,Lk=maxx∈xk,xk+1f,Nk=maxx∈xk,xk+1,x,u∈D∂f∂u,we have the following relation:(20)fξk,uξk−fck,uck≤ξk−ckMk+Nk⋅Lk.
Under the consideration ofξk,ck∈xk,xk+1, the following result is determined:(21)Qxk,xk+1,ξk,ck≤xk+1−xk2Mk+Nk⋅Lk.
The proof is completed.As shown in Theorem 1, it is found that the values ofMk, Nk, and Lk are dependent on the interval xk,xk+1. That is, a new function is constructed as follows:(22)Q¯kxk,xk+1=xk+1−xk2Mk+Nk⋅Lk.In addition, considering the error estimate of the approximate solutionuk+1, we give the following result.Theorem 2.
Assume that the functionfx,u is defined on D⊂R2 with fx,u:D⟶R. If fx,u is continuously differentiated on D, the error estimate is obtained as follows:(23)uxk+1−uk+1≤∑j=0kQ¯jxj,xj+1.Proof.
According to (14) and Theorem 1, we have(24)uxk+1−uk+1=uxk−uk+xk+1−xkfξk,uξk−fck,uck≤uxk−uk+Q¯kxk,xk+1≤uxk−1−uk−1+Q¯kxk,xk+1+Q¯k−1xk−1,xk⋯⋯⋯⋯⋯≤ux0−u0+Q¯kxk,xk+1+Q¯k−1xk−1,xk+⋯+Q¯0x0,x1=∑j=0kQ¯jxj,xj+1.
The proof is completed.Making use of Theorem 2, we can form a new optimization problem as follows:(25)minx0<x1<⋯<xk+1∑j=0kQ¯jxj,xj+1.Once the optimization problem (25) is solved, the nodes xjj=1,2,…,k+1 can be determined, and they are further used to determine uj through an Euler-type method such as one of (2)–(7). Hence, the approximate solution of ux x0≤x≤xk+1 is given as(26)uax=∑j=0k+1ujljx,where(27)l0x=x−x1x0−x1,x∈x0,x1,0,x∈x1,xk+1,ljx=x−xj−1xj−xj−1,x∈xj−1,xj,x−xj+1xj−xj+1,x∈xj,xj+1,0,x∉xj−1,xj+1,lk+1x=x−xkxk+1−xk,x∈xk,xk+1,0,x∈x0,xk.Of much interest is that the approximate solutionuk+1 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 uk+1 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 uk+1. 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 x0 to xk+1. The following Volterra integral equation is obtained:(28)uxk+1=u0+∫x0xk+1fx,udx.By considering all the grid pointsx0<x1<x2<⋯<xk+1, integral equation (28) is rewritten as(29)uxk+1=u0+∑j=0k∫xjxj+1fx,udx.Under the assumption of a continuous functionfx,u on a domain D, there is at least one ξj∈xj,xj+1 such that(30)uxk+1=u0+∑j=0kxj+1−xjfξj,uξj,where the integral mean value theorem has been applied [29]. When the values of ξj and uξj can be given explicitly, formula (30) is also explicit. However, it is always difficult to give the explicit values of ξj and uξj. Then, it is assumed that ξj and uξj are replaced by cj and ucj, respectively. Hence, the optimization problem is constructed as follows:(31)minuxk+1−uk+1=mincj∈xj,xj+1∑j=0kxj+1−xjfξj,uξj−fcj,ucj=mincj∈xj,xj+1∑j=0k∫xjxj+1fx,udx−∑j=0kxj+1−xjfcj,ucj=mincj∈xj,xj+1∑j=0k∫xjxj+1fx,udx−∑j=0k∫xjxj+1fcj,ucjdx=mincj∈xj,xj+1∑j=0k∫xjxj+1fx,u−fcj,ucjdx.Moreover, we introduce the Lipschitz conditions as follows:(32)fx,u−fcj,u≤L¯jx−cj,x,cj∈xj,xj+1,∀u,fx,ux−fx,ucj≤L˜jux−ucj,∀x∈xj,xj+1,where L¯j and L˜j are positive Lipschitz constants. It gives (33)∑j=0k∫xjxj+1fx,u−fcj,ucjdx=∑j=0k∫xjxj+1fx,u−fcj,ux+fcj,ux−fcj,ucjdx≤∑j=0k∫xjxj+1fx,u−fcj,ux+fcj,ux−fcj,ucjdx=∑j=0k∫xjxj+1L¯jx−cj+L˜jux−ucjdx.In addition, applying the differential mean value theorem, there is at least oneηj such that(34)ux−ucj=fηj,uηj⋅x−cj,where ηj is located between x and cj. By considering(35)Lj=maxx∈xj,xj+1f,we have(36)∑j=0k∫xjxj+1L¯jx−cj+L˜jux−ucjdx=∑j=0k∫xjxj+1L¯jx−cj+L˜jfηj,uηj⋅x−cjdx≤∑j=0k∫xjxj+1L¯j+Lj⋅L˜jx−cjdx.Applying the Cauchy–Schwarz inequality to (34) together with (33), the following relation is given:(37)∑j=0k∫xjxj+1fx,u−fcj,ucjdx≤∑j=0kL¯j+Lj⋅L˜j∫xjxj+1x−cjdx≤∑j=0kL¯j+Lj⋅L˜j∫xjxj+112dx⋅∫xjxj+1x−cj2dx1/2=∑j=0kL¯j+Lj⋅L˜jxj+1−xj3xj+1−cj3+xj−cj31/2.Then, the optimization problem (31) is transformed to the following form:(38)mincj∈xj,xj+1∑j=0kL¯j+Lj⋅L˜jxj+1−xj3xj+1−cj3+xj−cj31/2.It is easy to give the optimal solution of (38) as(39)Φ=36∑j=0kL¯j+Lj⋅L˜jxj+1−xj2,with cj=xj+xj+1/2. This means that when the optimization problem (38) is solved, the approximate solution of uxk+1 can be determined as follows:(40)uk+1=u0+∑j=0kxj+1−xjfxj+xj+12,uxj+xj+12.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:(41)F=minx0<x1<⋯<xk+1Φ.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 uk+1. 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 uxj+xj+1/2. For the sake of simplicity, the Euler formula is applied to give(42)uxj+xj+12≈uk+xj+1−xj2fxk,uk.Formula (40) is rewritten as(43)uk+1=u0+∑j=0kxj+1−xjfxj+xj+12,uj+xj+1−xj2fxj,uj.As compared with the differential method, the obtained functions (25) and (39) are similar, when the Lipschitz constants L¯j and L˜j are replaced by Mk and Nk, respectively.According to the above analysis, the two optimization problems (25) and (41) are constructed to determine the nonequivalent grids x0<x1<⋯<xk such that the approximate solution uk+1 to uxk+1 has more accuracy than that by using the equivalent grids. Furthermore, when we consider the grids a=x0<x1<⋯<xn=b and the error estimate on all the grids xk k=1,2,…,n, the optimization problem (13) can be reread as(44)min∑k=0n−1uxk+1−uk+1.In terms of (44) and (25), the derived optimization problem is obtained as(45)minx0<x1<⋯<xn∑k=0n−1∑j=0kQ¯jxj,xj+1.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]:(46)duxdx=x3−ux,x>1,u1=25,with the exact solution(47)ux=15x4+15x.
First, let us construct an optimization problem according to (25). In terms of (46), it gives(48)fx,u=x3−ux,meaning that(49)∂f∂x=3x2+ux2,∂f∂u=1x.
We further have the following estimates:(50)Mk=maxx∈xk,xk+13x2+ux2≤3xk+12+1xk2⋅maxx∈xk,xk+1u,Nk=maxx∈xk,xk+11x=1xk,Lk=maxx∈xk,xk+1x3−ux≤xk+13+1xkmaxx∈xk,xk+1u.
Then, it follows(51)Mk+Nk⋅Lk≤3xk+12+xk+13xk+2xk2maxx∈xk,xk+1u.
One can find from (51) that the value of maxx∈xk,xk+1u is unknown. Hence, the main term of (51) related to the nodes is extracted to yield the following function:(52)F=∑j=0kxj+1−xj23xj+12+xj+13xj+2xj2.
The optimization problem(53)minx0<x1<⋯<xk+1Fis used to determine the values of the nodes. Hence, the application of (43) leads to the computation formula as follows:(54)uk+1=u0+∑j=0kxj+1−xjxj+xj+123−2xj+1+xj3xj−xj+12xjuj+xj+1−xj2xj3.
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]:(55)v˗i+1=ωv˗i+U˗0,ϕ1⊗p˗i−x˗i+U˗0,ϕ2⊗p˗g−x˗i,x˗i+1=x˗i+v˗i.
The meanings of the terms in the above formulae are explained in Table1.
Moreover, for numerical computations, it is assumed thatω is monotonically decreasing from 0.9 to 0.4 with respect to the generation number and ϕi=2 i=1,2. When applying the PSO algorithm to the optimization problem (53), a particle is considered as the following k+2-dimensional vector:(56)x˗i=x0,x1,…,xk+1,subject to the constraint condition:(57)1=x0<x1<⋯<xk<xk+1=c,where c is a prescribed constant. The size of particle swarm is 1000, and the maximum of the generation is 200. The symbols p˗i and p˗g stand for the individual best position of a particle and the global best position of 1000 particles, respectively.
Figure1 is drawn to show the variation of the values of F versus the generation number with k=4 and xk+1=2. 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:(58)1.0000,1.2469,1.4630,1.6573,1.8352,2.0000.
The approximate solution using the nonequidistant grid points and formula (54) is obtained as u5=3.3046. When the equidistant grid points and (2) are used, we obtain u˜5=2.8407. The exact solution is u2=3.3. Clearly, the solution u5=3.3046 is more approximate to u2=3.3 than u˜5=2.8407. 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 k=4, we also choose xk+1=1.5,2.5,3.0,3.5 for computations, and the observations are shown in Table 2. It is found from Table 2 that the approximate solutions of ux5 determined by the Euler method are away from the exact one when the values of x5 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 x5. 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:(59)minx0<x1<⋯<xnFg,with(60)Fg=∑k=0n−1∑j=0kxj+1−xj23xj+12+xj+13xj+2xj2.
For the sake of comparisons, we choosen=5 and x5=2.0 for numerical computations. After running the PSO algorithm, the optimal solution is obtained, and the nonequivalent grids are given as(61)1.0000,1.1365,1.2849,1.4555,1.6703,2.0000.
The obtained approximation solution foru2 is u5=3.3068. Expectedly, it has more errors than the value u5=3.3046 in Table 2 by solving the optimization problem (53). For convenience, Figure 2 is drawn to show the variations of the absolute errors ux−ua1 and ux−ua2 versus the discrete grids, where ua1 and ua2 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 xn. The aim of the optimization problem (54) is the minimization of the global error on all discrete points.Table 1
The meanings of the terms in the formulae of particle swarm optimization.
x
→
i
The current position of particle
v
→
i
The current velocity of particle
p
→
i
The individual best position of particle
p
→
g
The global best position of all particles
U
→
0
,
ϕ
i
i
=
1,2
The vectors of random numbers uniformly distributed in0,ϕi
ω
The inertia weightFigure 1
The variations of the values ofF versus the generation number for k=4 and xk+1=2.Table 2
Comparisons with the typical Euler method.
The approximate solutionu5
x
5
Present study
Euler method
Exact solution
1.5
1.1465
1.0597
1.1458
2.0
3.3046
2.8407
3.3
2.5
7.9066
6.5183
7.8925
3.0
16.2985
13.0865
16.2667
3.5
30.1304
23.7542
30.0696Figure 2
The variations of the absolute errors versus the discrete points.Example 2.
In engineering applications, the widely studied Riccati differential equation is expressed as the following form [32, 33]:(62)duxdx=Pxu+Qxu2+Rx.
Here, we assume that it is subject to the initial value conditionu0=0 with Px=2, Qx=−1, and Rx=1. Then, the explicit solution is given as [34, 35](63)ux=1+2tanh2x+ln2−1.
Under the consideration of (46), we have(64)fx,u=2u−u2+1,∂f∂x=0,∂f∂u=2−2u.
In virtue of (25), an optimization problem is constructed as follows:(65)H=minx0<x1<⋯<xk+1∑j=0kxj+1−xj2.
For numerical computations, we choosex20=2 with k=19. 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 ui i=1,2,…,20 are obtained as(66)uk+1=u0+∑j=0kxj+1−xj1+2uj+xj+1−xj21+2uj−uj2−uj+xj+1−xj21+2uj−uj22.
In order to show the advantages of the proposed approach, the approximate solutions in [34, 35] are recalled and given in the following:(67)uex=x+x2+13x3−23x4+215x5,utx=10xx3−9x2+18x−602x5−31x4+200x3−630x2+720x−600.
As shown in Table3, the absolute errors uxk−uk, uxk−uexk, and uxk−utxk 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.Table 3
The optimized points and the absolute errors.
The present results
[34]
[35]
x
k
u
x
k
u
k
u
x
k
−
u
k
u
x
k
−
u
e
x
k
u
x
k
−
u
t
x
k
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.1149
0.1285
0.1277
0.0008
0.0000
0.0016
0.2316
0.2881
0.2865
0.0017
0.0005
0.0075
0.3408
0.4634
0.4612
0.0022
0.0016
0.0185
0.4552
0.6701
0.6677
0.0024
0.0023
0.0370
0.5660
0.8859
0.8837
0.0022
0.0003
0.0627
0.6790
1.1114
1.1097
0.0017
0.0106
0.0963
0.7910
1.3294
1.3282
0.0012
0.0326
0.1354
0.8989
1.5249
1.5239
0.0010
0.0670
0.1761
1.0072
1.7004
1.6993
0.0011
0.1140
0.2173
1.1091
1.8431
1.8417
0.0014
0.1659
0.2536
1.2136
1.9660
1.9643
0.0018
0.2210
0.2866
1.3113
2.0607
2.0587
0.0020
0.2675
0.3122
1.4064
2.1359
2.1338
0.0021
0.3012
0.3316
1.4967
2.1938
2.1916
0.0022
0.3167
0.3446
1.5820
2.2381
2.2360
0.0021
0.3119
0.3521
1.6680
2.2742
2.2723
0.0020
0.2840
0.3551
1.7498
2.3020
2.3002
0.0018
0.2329
0.3537
1.8369
2.3257
2.3241
0.0016
0.1497
0.3481
1.9208
2.3439
2.3425
0.0015
0.0400
0.3387
2.0000
2.3578
2.3565
0.0013
0.0911
0.3265


### 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:(68)Ld2Idt2+RdIdt+IC=0,t≥0,I0=0,dI0dt=10,where I is the electric current; L, R, and C are constants with R/L=20 and LC=100, respectively. The explicit solution can be computed as(69)It=25011113333e−100+31111/50t−e−100−31111/50t.
In order to apply the proposed method, it is assumed that(70)ut=It,dItdt=vt.
The initial problem (68) is rewritten as the system of the first order differential equations: (71)dutdt=vt,dvtdt=−utLC−RLvt,subject to the initial value conditions:(72)u0=0,v0=10.
In addition, we can also reexpress (71) as the style of the vector. That is, it is supposed that(73)Ut=utvt,dUtdt=u′tv′t.
Then, one has(74)dUtdt=AUt,t≥0,U0=U0,where(75)A=01−1LC−RL,U0=010.
For a series of grid points0=t0<t1<t2<⋯<tk+1, the application of (43) leads to(76)Uk+1=U0+∑j=0ktj+1−tjAUj+tj+1−tj2AUj,where Uj denotes the approximate solution of Utj. When the optimal solution Uk+1 is focused on, the grid points tj j=0,1,…,k+1 are determined by solving the following optimization problem:(77)G=mint0<t1<⋯<tk+1∑j=0ktj+1−tj2.
When considering the global optimal solution, the optimization problem for determining the grid pointstj j=0,1,…,n is constructed as follows:(78)Gg=mint0<t1<⋯<tn∑k=0n−1∑j=0ktj+1−tj2.
In what follows, the grid points0=t0<t1<t2<⋯<t5=0.1 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 as(79)0.0000,0.0200,0.0400,0.0600,0.0800,0.1000,0.0000,0.0088,0.0197,0.0343,0.0562,0.1000,respectively. 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 t. 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.Figure 3
The variations of the explicit and the approximate solutions. The approximate solutions with the marksE, G, and Gg mean that they are obtained by using the typical Euler method, the proposed method under solving the optimization problem G and that with solving the optimization problem Gg, respectively.


## 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]:(46)duxdx=x3−ux,x>1,u1=25,with the exact solution(47)ux=15x4+15x.
First, let us construct an optimization problem according to (25). In terms of (46), it gives(48)fx,u=x3−ux,meaning that(49)∂f∂x=3x2+ux2,∂f∂u=1x.
We further have the following estimates:(50)Mk=maxx∈xk,xk+13x2+ux2≤3xk+12+1xk2⋅maxx∈xk,xk+1u,Nk=maxx∈xk,xk+11x=1xk,Lk=maxx∈xk,xk+1x3−ux≤xk+13+1xkmaxx∈xk,xk+1u.
Then, it follows(51)Mk+Nk⋅Lk≤3xk+12+xk+13xk+2xk2maxx∈xk,xk+1u.
One can find from (51) that the value of maxx∈xk,xk+1u is unknown. Hence, the main term of (51) related to the nodes is extracted to yield the following function:(52)F=∑j=0kxj+1−xj23xj+12+xj+13xj+2xj2.
The optimization problem(53)minx0<x1<⋯<xk+1Fis used to determine the values of the nodes. Hence, the application of (43) leads to the computation formula as follows:(54)uk+1=u0+∑j=0kxj+1−xjxj+xj+123−2xj+1+xj3xj−xj+12xjuj+xj+1−xj2xj3.
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]:(55)v˗i+1=ωv˗i+U˗0,ϕ1⊗p˗i−x˗i+U˗0,ϕ2⊗p˗g−x˗i,x˗i+1=x˗i+v˗i.
The meanings of the terms in the above formulae are explained in Table1.
Moreover, for numerical computations, it is assumed thatω is monotonically decreasing from 0.9 to 0.4 with respect to the generation number and ϕi=2 i=1,2. When applying the PSO algorithm to the optimization problem (53), a particle is considered as the following k+2-dimensional vector:(56)x˗i=x0,x1,…,xk+1,subject to the constraint condition:(57)1=x0<x1<⋯<xk<xk+1=c,where c is a prescribed constant. The size of particle swarm is 1000, and the maximum of the generation is 200. The symbols p˗i and p˗g stand for the individual best position of a particle and the global best position of 1000 particles, respectively.
Figure1 is drawn to show the variation of the values of F versus the generation number with k=4 and xk+1=2. 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:(58)1.0000,1.2469,1.4630,1.6573,1.8352,2.0000.
The approximate solution using the nonequidistant grid points and formula (54) is obtained as u5=3.3046. When the equidistant grid points and (2) are used, we obtain u˜5=2.8407. The exact solution is u2=3.3. Clearly, the solution u5=3.3046 is more approximate to u2=3.3 than u˜5=2.8407. 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 k=4, we also choose xk+1=1.5,2.5,3.0,3.5 for computations, and the observations are shown in Table 2. It is found from Table 2 that the approximate solutions of ux5 determined by the Euler method are away from the exact one when the values of x5 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 x5. 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:(59)minx0<x1<⋯<xnFg,with(60)Fg=∑k=0n−1∑j=0kxj+1−xj23xj+12+xj+13xj+2xj2.
For the sake of comparisons, we choosen=5 and x5=2.0 for numerical computations. After running the PSO algorithm, the optimal solution is obtained, and the nonequivalent grids are given as(61)1.0000,1.1365,1.2849,1.4555,1.6703,2.0000.
The obtained approximation solution foru2 is u5=3.3068. Expectedly, it has more errors than the value u5=3.3046 in Table 2 by solving the optimization problem (53). For convenience, Figure 2 is drawn to show the variations of the absolute errors ux−ua1 and ux−ua2 versus the discrete grids, where ua1 and ua2 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 xn. The aim of the optimization problem (54) is the minimization of the global error on all discrete points.Table 1
The meanings of the terms in the formulae of particle swarm optimization.
x
→
i
The current position of particle
v
→
i
The current velocity of particle
p
→
i
The individual best position of particle
p
→
g
The global best position of all particles
U
→
0
,
ϕ
i
i
=
1,2
The vectors of random numbers uniformly distributed in0,ϕi
ω
The inertia weightFigure 1
The variations of the values ofF versus the generation number for k=4 and xk+1=2.Table 2
Comparisons with the typical Euler method.
The approximate solutionu5
x
5
Present study
Euler method
Exact solution
1.5
1.1465
1.0597
1.1458
2.0
3.3046
2.8407
3.3
2.5
7.9066
6.5183
7.8925
3.0
16.2985
13.0865
16.2667
3.5
30.1304
23.7542
30.0696Figure 2
The variations of the absolute errors versus the discrete points.Example 2.
In engineering applications, the widely studied Riccati differential equation is expressed as the following form [32, 33]:(62)duxdx=Pxu+Qxu2+Rx.
Here, we assume that it is subject to the initial value conditionu0=0 with Px=2, Qx=−1, and Rx=1. Then, the explicit solution is given as [34, 35](63)ux=1+2tanh2x+ln2−1.
Under the consideration of (46), we have(64)fx,u=2u−u2+1,∂f∂x=0,∂f∂u=2−2u.
In virtue of (25), an optimization problem is constructed as follows:(65)H=minx0<x1<⋯<xk+1∑j=0kxj+1−xj2.
For numerical computations, we choosex20=2 with k=19. 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 ui i=1,2,…,20 are obtained as(66)uk+1=u0+∑j=0kxj+1−xj1+2uj+xj+1−xj21+2uj−uj2−uj+xj+1−xj21+2uj−uj22.
In order to show the advantages of the proposed approach, the approximate solutions in [34, 35] are recalled and given in the following:(67)uex=x+x2+13x3−23x4+215x5,utx=10xx3−9x2+18x−602x5−31x4+200x3−630x2+720x−600.
As shown in Table3, the absolute errors uxk−uk, uxk−uexk, and uxk−utxk 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.Table 3
The optimized points and the absolute errors.
The present results
[34]
[35]
x
k
u
x
k
u
k
u
x
k
−
u
k
u
x
k
−
u
e
x
k
u
x
k
−
u
t
x
k
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.1149
0.1285
0.1277
0.0008
0.0000
0.0016
0.2316
0.2881
0.2865
0.0017
0.0005
0.0075
0.3408
0.4634
0.4612
0.0022
0.0016
0.0185
0.4552
0.6701
0.6677
0.0024
0.0023
0.0370
0.5660
0.8859
0.8837
0.0022
0.0003
0.0627
0.6790
1.1114
1.1097
0.0017
0.0106
0.0963
0.7910
1.3294
1.3282
0.0012
0.0326
0.1354
0.8989
1.5249
1.5239
0.0010
0.0670
0.1761
1.0072
1.7004
1.6993
0.0011
0.1140
0.2173
1.1091
1.8431
1.8417
0.0014
0.1659
0.2536
1.2136
1.9660
1.9643
0.0018
0.2210
0.2866
1.3113
2.0607
2.0587
0.0020
0.2675
0.3122
1.4064
2.1359
2.1338
0.0021
0.3012
0.3316
1.4967
2.1938
2.1916
0.0022
0.3167
0.3446
1.5820
2.2381
2.2360
0.0021
0.3119
0.3521
1.6680
2.2742
2.2723
0.0020
0.2840
0.3551
1.7498
2.3020
2.3002
0.0018
0.2329
0.3537
1.8369
2.3257
2.3241
0.0016
0.1497
0.3481
1.9208
2.3439
2.3425
0.0015
0.0400
0.3387
2.0000
2.3578
2.3565
0.0013
0.0911
0.3265


## 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:(68)Ld2Idt2+RdIdt+IC=0,t≥0,I0=0,dI0dt=10,where I is the electric current; L, R, and C are constants with R/L=20 and LC=100, respectively. The explicit solution can be computed as(69)It=25011113333e−100+31111/50t−e−100−31111/50t.
In order to apply the proposed method, it is assumed that(70)ut=It,dItdt=vt.
The initial problem (68) is rewritten as the system of the first order differential equations: (71)dutdt=vt,dvtdt=−utLC−RLvt,subject to the initial value conditions:(72)u0=0,v0=10.
In addition, we can also reexpress (71) as the style of the vector. That is, it is supposed that(73)Ut=utvt,dUtdt=u′tv′t.
Then, one has(74)dUtdt=AUt,t≥0,U0=U0,where(75)A=01−1LC−RL,U0=010.
For a series of grid points0=t0<t1<t2<⋯<tk+1, the application of (43) leads to(76)Uk+1=U0+∑j=0ktj+1−tjAUj+tj+1−tj2AUj,where Uj denotes the approximate solution of Utj. When the optimal solution Uk+1 is focused on, the grid points tj j=0,1,…,k+1 are determined by solving the following optimization problem:(77)G=mint0<t1<⋯<tk+1∑j=0ktj+1−tj2.
When considering the global optimal solution, the optimization problem for determining the grid pointstj j=0,1,…,n is constructed as follows:(78)Gg=mint0<t1<⋯<tn∑k=0n−1∑j=0ktj+1−tj2.
In what follows, the grid points0=t0<t1<t2<⋯<t5=0.1 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 as(79)0.0000,0.0200,0.0400,0.0600,0.0800,0.1000,0.0000,0.0088,0.0197,0.0343,0.0562,0.1000,respectively. 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 t. 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.Figure 3
The variations of the explicit and the approximate solutions. The approximate solutions with the marksE, G, and Gg mean that they are obtained by using the typical Euler method, the proposed method under solving the optimization problem G and that with solving the optimization problem Gg, respectively.


## 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 degreeIt 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.



---

*Source: 9071236-2019-12-12.xml*