The Integrating Factors for Riccati and Abel Di ff erential Equations

We can recast the Riccati and Abel differential equations into new forms in terms of introduced integrating factors. Therefore, the Lie-type systems endowing with transformation Lie-groups S L(2,R) can be obtained. The solution of second-order linear homogeneous differential equation is an integrating factor of the corresponding Riccati differential equation. The numerical schemes which are developed to fulfil the Lie-group property have better accuracy and stability than other schemes. We demonstrate that upon applying the group-preserving scheme (GPS) to the logistic differential equation, it is not only qualitatively correct for all values of time stepsize h, and is also the most accurate one among all numerical schemes compared in this paper. The group-preserving schemes derived for the Riccati differential equation, second-order linear homogeneous and non-homogeneous differential equations, the Abel differential equation and higher-order nonlinear differential equations all have accuracy better than O(h2).


Introduction
The Lie-groups play an important role to understand the geometric structure of differential equations, which are very useful in the development of some superior numerical methods to integrate ordinary differential equations (ODEs), and to retain the invariant property of dynamical system.By sharing the geometric structure and invariant with the original ODEs, the new methods are more accurate, more stable and more effective than conventional numerical methods.In the last few decades there has been a substantial development in the geometric integrators to solve ODEs evolving on the Lie groups (Munthe-Kass, 1999;Iserles et al., 2000;Liu, 2007;Hochbruck & Ostermann, 2010).The Lie-group schemes find numerical approximations to the solution of where the exact solution Y is a matrix Lie group with A being a matrix function on the associated Lie algebra.
In this paper we introduce the integrating factor methods for the Riccati, Abel and higher-order nonlinear differential equations, and then use a Lie-group S L(2, R) approach for the solutions of the above equations.The Lie-group S L(2, R) has been used to develop the shooting method to solve the generalized Sturm-Liouville problem by Liu (2012), to solve a singular ϕ-Laplacian nonlinear ODE by Liu (2013a), and to solve a nonlinear heat transfer equation by Hashemi (2015).
For the first-order linear ODE: ẏ(t) + p(t)y(t) = q(t), (1) we have already learned that the integrating factor is defined by which can transform Equation (1), after multiplying both the sides by I(t), into a total differential form on the left-hand side: d dt [I(t)y(t)] = I(t)q(t). (3) Such that the general solution of Equation ( 1) can be obtained simply by integrating the above equation.
The key point by passing from Equation (1) to Equation ( 3) is the idea of integrating factor, leading to a total differential form in Equation ( 3).This idea can be extended to the Riccati equation.It not only provides a total differential form of the left-hand side and also linearizes the nonlinear differential equation exactly, of which the integrating factor u(t) is a solution of a second-order linear ODE: In the course of Engineering Mathematics we may teach the students to write Equation (4) as a system of two first-order ODEs by d dt However, letting be the integrating factor of Equation ( 4) we have It shows again that the left-hand side can be managed to a total differential form by using the concept of integrating factor.Furthermore, from Equation ( 7) we have a self-adjoint system: This system is better than system (5), because it allows an internal symmetry of S L(2, R).
Although, the author (Liu, 2001;Liu, 2013b;Liu, 2015) has previously developed several general Lie-group theories for the general nonlinear dynamical system, the specialization to the Lie-group S L(2, R) for the secondorder linear ODE and Riccati equation and finding that the solution of the former is an integrating factor of the latter is interesting, which deserves a further attention by sharing the Lie-group property to develop more efficient numerical scheme for these two equations.
There are many ODEs that the concept of integrating factor can be used to derive new forms of ODEs, of which the internal symmetry appears naturally.Liu (2001) was the first by using the technique of integrating factor to transform the nonlinear ODEs into an augmented system in the Minkowski space, where the length of state vector plays the role of integrating factor.Then, Liu (2006) introduced the concept of multiple integrating factors to solve the ODEs with multiple constraints.
Many numerical methods have been developed for the general purpose of integrating the system of ODEs to obtain numerical solutions.Among those numerical methods, the most known are the Euler method and the Runge-Kutta method, but less known is the group-preserving scheme (Liu, 2001).These methods discretize the differential equations system to produce difference equations or mappings.The different methods obtain different mappings from the same differential equations, but they have the same aim in that the dynamics of the mapping should correspond as closely as possible to the dynamics of original differential equations.However, it is a major question whether the corresponding approximation has the same behavior as the original differential equations.In recent years, a major area is to study the difference in the exact dynamics of differential equations and the dynamics of the corresponding discretized mappings.It is found that the most serious problems for the numerical discretizations are the appearance of so-called spurious solutions and ghost fixed points in the discretized dynamics.For example, Chen and Solis (1998) have investigated the appearance of spurious solutions when some first-order ODEs are discretized by using the Runge-Kutta methods, and showed that the resulting schemes may produce unrelated bifurcation phenomena.
Of all the first-order nonlinear differential equations, the Riccati differential equation occupies perhaps the most important place, because it is closely related to the second-order linear homogeneous differential equations, which appear in many physical applications.In the Engineering Mathematics we have learned some special functions and series solutions for the Bessel equation, the Legendre equation, the Airy equation, etc., but from the computational point of view it is not convenient to use these methods to calculate the solutions.In this paper we attempt to develop suitable numerical schemes for the Riccati differential equation, second-order linear differential equations, and the Abel differential equation.Although the above systems will be shown exhibiting certain symmetry, there has no development of the numerical scheme to preserve it.Really, there have many issues concerning the improvement of the quality of solution by producing numerical solutions of ODEs that possess symmetries and structures and also utilize the symmetries to ensure the computational stability.
The arrangements of this paper are given as follows.In Section 2 we introduce a Lie-group symmetry for the Riccati differential equation, and then we explain that its integrating factor is a solution of the second-order linear differential equation in Section 3. Based on the Lie-symmetry we develop a group-preserving scheme for the Riccati differential equation in Section 4, and the group-preserving scheme for the logistic differential equation is developed in Section 5, which is compared with other numerical schemes.The group-preserving scheme for secondorder linear differential equation is developed in Section 6, while that for second-order linear non-homogeneous differential equation is developed in Section 7. The group-preserving schemes for the Abel differential equation and other higher-order nonlinear differential equations are developed in Section 8. Finally, the conclusions are drawn in Section 9.

Symmetry Group
The following Riccati differential equation: can be written as Upon defining the integrating factor: Equation ( 10) becomes On the other hand, from Equation ( 11) it follows that Equations ( 12) and ( 13) together constitute a linear system: where We thus represent the Riccati differential equation in the projective space (uy, u) as a system of two first-order linear ODEs.
The above coefficient matrix A satisfying tr A = 0 is a Lie algebra of S L(2, R), where tr denotes the trace.Thus, the one-parameter subgroup generated by A gives the following transformation formula for (uy, u): where G is an element of the two-dimensional special linear group, satisfying det in which det stands for the determinant, and I 2 is a second order identity matrix.The above fact prompts us to develop a suitable numerical scheme to preserve the Lie-group property.

The Second-order Linear Differential Equations
It follows from Equations ( 13) and ( 12) a second-order linear ODE for u: where For a given a(t) it follows directly that On the other hand, without lost any generality we can simply take If G is a fundamental matrix of Equation ( 14) satisfying then from Equation ( 16), with its uy replaced by u/p due to q(t) = 0 in Equation ( 13), we can obtain Noting that G 21 (t) and G 22 (t) are linearly independent solutions of Equation ( 18), because in Equation ( 23), is a general solution of Equation ( 18).
Differentiating the second equation with respect to t and comparing it with the first equation, we obtain Substituting the above two equations for G 11 (t) and which is known as the Abel identity.It is clear that such an identity is a direct result of the Lie-group property (17).

Numerical Scheme for Riccati Differential Equation
A numerical scheme based on the Lie-group property can be utilized to enhance the computational accuracy and efficiency.The time-centered Euler scheme for Equation ( 14) is where Here n ∈ Z + , u n and y n denote respectively the value of u and y at the discrete time t n , and τ is one half of the stepsize h, that is, τ := h/2 = ∆t/2.The notation n is used to denote n + 1/2, and A(n) means that A takes value at t n + τ.From Equation (15) we have Substituting the above A(n) into Equation ( 27) and through some manipulations we find that A straightforward calculation shows that det G(n) = 1 holds.Thus the above numerical scheme meets the requirement for preserving the Lie-group property.Substituting Equation ( 29) for G(n) into Equation ( 26) and taking the projection, we obtain which is a group-preserving scheme (GPS) for the Riccati differential equation ( 9).
For the purpose of comparison we give a numerical example to calculate the solution of Equation ( 9) with p(t) = −2 cos t, q(t) = 1/2+2t cos t, and r(t) = 1+t+2t 2 cos t, of which the exact solution is y(t) = t−[sin t−cos t+c exp t] −1 , where c = 1 − 1/y(0).Starting from the initial condition y(0) = −1, using h = 0.001 and substituting the above coefficient functions into Equation ( 30) we can obtain the numerical solution, which upon comparing with the closed-form solution we show the solution error, which is obtained by substracting the numerical solution from the closed-form solution, in Figure 1.It can be seen that the accuracy of scheme ( 30) is better than the order of O(h 2 ).
From Equation ( 30) we can conclude that a suitable scheme for the Riccati differential equation is a linear fraction map of the following form:

Numerical Schemes for Logistic Differential Equation
Now, let us apply scheme (30) to the logistic differential equation: whose solution is given by y with y(0) = y 0 .Inserting p = λ, q = −λ/2 and r = 0 into Equation ( 30), we obtain a discretized map for the logistic differential equation: Obviously, g has only two fixed points y = 0 and y = 1, and because of the first fixed point is unstable and the second one is stable when λ > 0. Conversely, lead to the fixed point y = 0 to be stable and the fixed point y = 1 to be unstable when λ < 0. Thus, the map (35) has the same properties as the original equation ( 33), no matter what λ is.Furthermore, by the induction method the solution of Equation ( 35) is found to be It has the same form as Equation ( 34), and no matter what the initial value y 0 is and what λ is, it tends to the stable fixed point.
In Figure 2, where λ = 2, h = 0.001 and y 0 = 0.1 were used, we compare the numerical solution with the closedform solution in terms of solution error.The accuracy of scheme ( 35) can be seen almost in the order of O(h 3 ).Also the comparisons with the following schemes are made in Figure 2: Here, the first two schemes are proposed by Wang et al. (1992), and the third is the combination of the first two schemes (Gander & Meyer-Spasche, 2000, p.224).The fourth is called by Mickens (1999) the nonstandard finite-difference scheme, the factor 1 − exp(−λh) in which is known as the denominator function.
Th accuracies of the above numerical schemes are compared in Figure 2, which indicate that scheme (39) has accuracy better than the order of O(h 2 ), the accuracy of shemes (37) and ( 38) is better than the order of O(h), but scheme (40) has accuracy in the order of O(h).As shown in Figure 2, schemes (37) and (38) lead to the same absolute solution error but with difference in sign, and so we introduce a scheme being the average of them: whose accuracy is better than schemes (37) and ( 38), but is worse than scheme (39).Now we discuss something about the group-preserving scheme.When λh < 1, scheme (37) can be written as the form of Equation ( 31) with It satisfies condition (32), so that scheme (37) is a group-preserving scheme in the range λh < 1.However, in the range of λh > 1, scheme (37) fails to be a group-preserving scheme.Gander & Meyer-Spasche (2000) have proved that scheme (37) behaves qualitatively correctly for λh < 1.For 1 < λh < 2, it converges to the correct limit oscillatorily.For λh > 2, it behaves wrongly.Thus, a scheme which can preserve the Lie-group property is very important.Scheme ( 38) is a group-preserving scheme if λh > −1; otherwise, it fails to be a group-preserving scheme.Similarly, scheme (39) is a group-preserving scheme in the range λh < 2, and when λh > 2, it fails to be a group-preserving scheme.It is obvious that only the scheme ( 35) is a group-preserving for all λ.Schemes ( 40) and ( 41) are not of the linear fraction type as Equation ( 31) is, and thus are not of the group-preserving type.
There have several criteria to assess the performance of a numerical scheme, among, the order of accuracy, stability property and preserving the dynamics of the original ODEs are the most important three factors.In this regard, the above introduced group-preserving scheme ( 35) is the best one for the discretizations of the logistic differential equation (33).

Numerical Scheme for Second-order Linear Differential Equations
In order to find the numerical scheme for Equation ( 18), we let u = p(t)uy by Equation ( 13) due to q(t) = 0; hence, from Equations ( 26) and ( 29) the numerical scheme for u and u can be written as where p is defined in Equation ( 21).The above scheme is suitable to calculating the solution of the second-order linear differential equation ( 18), whose initial values of u and u are specified.
Let us consider a special case ü + b(t)u = 0. ( 43) The general second-order linear homogeneous differential equations can always be transformed into the above type.For this case let p = 1 in Equation ( 42) we obtain Through some manipulations we obtain It bears certain similarity with the following two schemes: developed by Mickens & Ramadhani (1992), and generalized from the Mickens-Ramadhani scheme by Chen et al. (1993).
In order to compare the above three numerical schemes we let b(t) = 1/t 2 in the range t ∈ [1, 10].For such b(t), Equation ( 43) has closed-form solution as follows: The solution error of the above three schemes are compared in Figure 3, where u(1) = u(1) = 1 were taken, and the stepsize h = 0.001 was used.As denoted by Chen et al. (1993) the modified scheme ( 47) is of O(h 4 ) and the Mickens-Ramadhani scheme ( 45) is of O(h 2 ).It can be seen that the accuracy of scheme ( 45) is between schemes ( 46) and ( 47).However the latter two schemes do not share the Lie-group property when b(t) is a time-varying function as this case is.
When we consider b(t) = ω 2 for the linear harmonic oscillator equation, the numerical scheme in Equation ( 44) is simplified to Here we require τ 2 ω 2 < 1.It deserves to note that the above scheme conserves the total energy On the other hand, when applying the nonstand difference scheme of Mickens to Equation ( 49) gives Refer Equations (2.40a) and (2.40b) in (Mickens, 2000).Obviously, Equation ( 51) also holds for this scheme.
Comparing the coefficients matrices in schemes ( 50) and ( 52), we note that both the transition matrices are of S L(2, R), and make the above two schemes slightly different in the order of O(h 4 ).

Numerical Scheme for Second-order Linear Non-homogeneous Differential Equations
Now, let us consider the following second-order linear non-homogeneous differential equation: with a forcing term f (t).It is equivalent to solve the following non-homogeneous linear system: The solution is given by where u(t i ) and u(t i ) are the initial values prescribed at time t i .The second line in Equation ( 55) is really a particular solution of Equation ( 53).It leads to the same result as that obtained by the method of variation of parameters.In Equation ( 54) the first equation d dt is known as the non-homogeneous Sturm-Liouville equation, and the differential operator is self-adjoint.The usual way to systemize Equation ( 53) is given by However, the system matrix of this equation does not has any special structure in it.Thus, the numerical method based on this equation will not capture the main Lie-group property inherent in such an equation.
For the non-homogeneous equation ( 53) we can let t → t n+1 , t i → t n , and apply the trapezoidal quadrature to the integrals in Equations ( 55) and ( 56) to obtain where In order to investigate the accuracy of the above numerical scheme we consider the following differential equation: in the range t ∈ [1, 10], with u(1) = 2 and u(1) = 1.Equation ( 60) has a closed-form solution given as follows: where Substituting p(t) = 1, b(t) = 1/t 2 and f (t) = ln t/t 2 into the above scheme we can calculate u, whose error is shown in Figure 4, where the stepsize h = 0.001 was used in the calculation.It can be seen that the accuracy is better than the order of O(h 2 ).
8. Abel Differential Equation and Higher-order Nonlinear ODE

The Integrating Factor for Abel Differential Equation
A direct extension of the Riccati differential equation is the Abel differential equation (Sachdev, 1991): which can be written as dy dt + [q(t) + p(t)y + e(t)y 2 ]y = r(t) − q(t)y. (63) Upon defining the integrating factor: Equation ( 62) becomes On the other hand, from Equation ( 64) it follows that du dt = q(t)u + p(t)uy + e(t)uy 2 . (66) Equations ( 65) and (66) together can be written as where Because the coefficient matrix A satisfies trA = 0, the resulting Lie-group is S L(2, R), which is a representation of the Abel differential equation in the projective space (uy, u) as a system of two first-order differential equations.However, the above A is different from that in Equation ( 15) with an extra term e(t)y = e(t)uy/u, which renders Equation (67) a quasi-linear system endowing with a local Lie-algebra element of sl(2, R).

The S L(2, R) Group-preserving Scheme
According to Equation (67), we can develop an implicit scheme based on S L(2, R) for the integration of Equation ( 62), of which one can view Āk := [ − qk rk pk + ēk ȳk qk ] as a constant matrix within a small time increment, where tk = (1 − θ)t k + θt k+1 , qk = q( tk ), rk = r( tk ), pk = p( tk ), ēk = e( tk ) and ȳk = (1 − θ)y k + θy k+1 .The corresponding Lie-group G k ∈ S L(2, R) is obtained by using the exponential mapping: Let For the special case with r(t) = 0 we have For the case r(t) 0, through some derivations we can obtain In above, This scheme is implicit because it also depends on y k+1 , which requires an iteration to determine the value of (z k+1 := u k+1 y k+1 , u k+1 ) at the next time step.The numerical process of group-preserving scheme (GPS) is given as follows.
For the Abel differential equation, there are only rare cases that the analytical solutions exist as summarized in Murphy (1960).In general the solutions are implicit.For example, dy dt = a 0 y + b 0 y 3 , y(0) = y 0 (77) has an analytical solution: where C = a 0 /y 2 0 + b 0 .We fix a 0 = 2, b 0 = 1 and y 0 = 1, and compare the solutions obtained by the GPS with the analytical solutions in Table 1.In order to test the accuracy of the above GPS we let y = sin t be an exact solution.Substituting p(t) = 0.5, q(t) = 1/2, e(t) = −1, y = sin t and ẏ = cos t into Equation (62) we can obtain r(t).Under h = 0.001 and ε = 10 −10 we use the above GPS to calculate y, whose number of iterations and the absolute error of numerical solution are shown in Figure 5.Although under a stringent convergence criterion ε = 10 −10 , the number of iterations is either three or four.It can be seen that the accuracy is better than the order of O(h 2 ).
Let us consider (Adeyeye et al., 2015) which is a special case of Equation ( 62) with p(t) = 0, q(t) = 1/2, e(t) = −1, and r(t) = 1.Under h = 0.001 and ε = 10 −10 we use the GPS to calculate y, whose number of iterations is shown in Figure 6(a), and the solution is compared with that calculated by the fourth-order Runge-Kutta method (RK4) in Figure 6 We must emphasize that the results shown in (Adeyeye et al., 2015;Jaradat, 2008) are incorrect after t = 1.2, whose errors are over 0.1 and more large up to 17.5 when time tends to t = 1.92,where the GPS leads to the value 21.22595.

The Extension to Highly Nonlinear ODEs
Through the above studies we can extend the Lie-group approach to highly nonlinear ODE: (b), and the difference is shown in Figure 6(c).Up to t = 1.92 the solution fast tends to large value 21.22595, which causes the numerical solution obtained by the RK4 blowing up as shown in Figure 6(b) by dashed line.However, the GPS is also available up to t = 1.92 as shown in Figure 6(b) by solid line.

Table 1 .
Comparing exact and GPS solutions for an Abel equation