Multi-point Taylor series to solve differential equations

The use of Taylor series is an effective numerical method to solve ordinary differential equations but this fails when the sought function is not analytic or when it has singularities close to the domain. These drawbacks can be partially removed by considering multi-point Taylor series, but up to now there are only few applications of the latter method in the literature and not for problems with very localized solutions. In this respect, a new numerical procedure is presented that works for an arbitrary cloud of expansion points and it is assessed from several numerical experiments.


1.
Introduction. It is quite common to solve Ordinary Differential Equations (ODE) from Taylor series [20,4] but much less for Partial Differential Equations (PDE). A new numerical technique based on Taylor series, named Taylor Meshless Method, has been presented rather recently [24,21,22,23] to solve elliptic PDE's: the unknowns and the residual are written in the form of high order polynomials. To account for the PDE, a truncated Taylor series of the residual is vanished, what greatly reduces the number of active degrees of freedom (dof). For instance a boundary value problem has been solved accurately with only 90 dof's while more than 40000 are needed to get the same result with linear triangular finite elements [21].
Nevertheless, as it is natural within Taylor series method, the convergence of the method depends strongly on the singularities of the unknown. Now singularities are present in practical engineering computations, because of domains with corners, pointwise forces, unilateral contact. . . Thus the treatment of such problems is a challenge and an interesting answer is offered by multi-domain approaches [22], each subdomain being associated with a Taylor series. Multi-point Taylor series could be an alternative way to circumvent these singularities and this is discussed in this paper in the case of ODE's.

DJÉDJÉ SYLVAIN ZÉZÉ, MICHEL POTIER-FERRY AND YANNICK TAMPANGO
The origin of multi-point Taylor approximations seems rather old [3] and this has been applied mainly to get accurate estimates of special functions like hypergeometric functions, see for instance [10,11,12]. The two-point Taylor expansions of an analytic function have been established in [13] where explicit forms of the coefficients are given. When two Taylor series at order p around points z 1 and z 2 are given, a new polynomial of degree 2p + 1 is deduced and this two-point Taylor expansion can converge even if there is a close singularity. The same holds with several arbitrarily located expansion points [14]. There are only few attempts of applications to boundary value problems for ODE's [15,5] and 2D second order partial differential equations [16]. In [15], 2-point and 3-point Taylor expansions provide accurate solutions in cases where the 1-point Taylor series cannot converge, according to theoretical results. So this technique seems promising in view of applications to boundary value problems of practical interest.
In this paper, a generic procedure is proposed to define multi-point Taylor series for any cloud of expansion points, by simply combining one-point Taylor series with the Hermite-Birkhoff interpolation that is well known in numerical analysis [18]. According to the framework of true meshless methods [9], the discretisation is defined only from points: collocation points and points to define shape functions satisfying the differential equation. Several numerical tests will be carefully discussed to assess the ability of multi-point Taylor series as compared to one-point series and or to multi-domain approaches. Meshless methods [17] were introduced to bypass the difficulty of re-meshing because it is easier to improve the approximation by adding points in well chosen regions. This implies to define robust error criteria to find places where an enrichment is needed [8]. Such an approach will be sketched in the case of multi-point Taylor series.

2.
Taylor series and two-point boundary value problems. Let us consider a two-point boundary value problem as where the given functions f(x) and g(x) are smooth. Well-posed Initial Value Problems can be associated to the differential equation and their solutions can be defined as functions of the two initial data u 0 , u 1 . As usual, the function R(u(x)) is called residual. If all these functions f, g, and u are written in the form of Taylor series, a simple recurrence yields all the derivatives of the unknown u(x) at x = x 1 as a function of u 0 , u 1 [20]. The idea is to combine several Taylor series with the help of the well known Hermite-Birkhoff interpolation [18] that defines a high order interpolation polynomial from all the derivatives of u(x) up to order "p" at a cloud of points (x 1 , x 2 , x 3 , . . . , x n ). These latter points can be called expansion points or high order collocation points. This leads to a family of polynomials depending of these 2n parameters. In the case of a one-point Taylor series (n = 1), these two parameters are determined by the two boundary conditions. If one applies a multi-point Taylor series (n > 1), the boundary conditions are not sufficient and one has to go back to the differential equation and to account for more information about the residual. In [15], these last equations follow from higher order derivatives of the residual at the same cloud of points. In this paper, we prefer to vanish the residual at a secondary cloud of points (y 1 , y 2 , y 3 , . . . , y 2n−2 ), which seems more flexible and avoids to compute high order derivatives of the interpolation polynomials. In summary, two methods will be used to discretize the differential equations : the derivatives of the residual up order "p" are put to zero at the main cloud of expansion points (x 1 , x 2 , x 3 , . . . , x n ) and only the residual at the secondary cloud of collocation points (y 1 , y 2 , y 3 , . . . , y 2n−2 ). If necessary, the problem can be solved in a piecewise manner with a multi-point series per subdomain. We are now describing the details of the algorithm.
2.1. Basic one-point Taylor series. The unknown u(x) is sought in the form of a truncated Taylor series of degree p: The input functions f (x) and g(x) and the residual R(x) are written in the same manner. There is an obvious connection between the derivatives and the coefficients of the series: j!u j = d j u dx j (x c ). Thus a Taylor series like (3) can be characterized by a vector bringing together all the coefficients of the series: Classically, when using Taylor series, the initial value problem (2) is solved by setting to zero the derivatives of the residual up to the order p-2: The latter equations (5) permit to solve the Initial Value Problem (2) because all the derivatives of u(x) can be obtained from u 0 and u 1 recursively. Hence the unknown vector U is an affine function of u 0 and u 1 and can be written in the form: The components of the vectors α, β, γ are easily deduced from the recurrence formula (5), α by starting from u 0 = 1, u 1 = 0, g(x) = 0, β by starting from u 0 = 0, u 1 = 1, g(x) = 0 , γ by starting from u 0 = 0, u 1 = 0, g(x). When this onepoint series is applied to the boundary value problem in the whole domain [−L, L], the two remaining unknowns u 0 and u 1 can be obtained from the two boundary conditions u(−L) = u(L) = 0.

2.2.
A general procedure to construct multi-point Taylor series. Hermite interpolation is very classical to build shape functions for beam bending problems: typically one defines polynomials of degree 3 from the values of the function u(x) and its derivatives u (x) at two nodes [19]. It is known that this simple interpolation can be extended to any degree p and to any cloud of points, see [18]. This cloud of n points will be denoted by x 1 , x 2 , · · · , x n . In each point of the cloud, one defines the vector of the Taylor coefficients as in (4): After solving numerically the equations (5) at the point x i , the vector U i is expressed as a function of the two variables u i0 , u i1 as in (6): where (α i ) j denotes the j th component of the vector α i (of course the same holds for the other vectors β i , γ i ).
Next one applies the Hermite-Birkoff interpolation method that defines a polynomial of degree n(p + 1) − 1 from the derivatives U i at the chosen cloud of points: The polynomials L ij (x) can be built by a classical procedure that is recalled in the Appendix A. One can combine the result of the Taylor method (8) with the interpolation formula (9), which permits to express the interpolated polynomial u(x) as a function of the values of the function and its first derivative on the cloud of points. For convenience, these 2n variables are aggregated in a vector that will be the unknown of the discrete problem: t q = (u 10 u 20 · · · u n0 u 11 u 21 · · · u n1 ) = (q 1 q 2 · · · q n q n+1 q n+2 · · · q 2n ) ∈ R 2n From (8), (9), one gets Hence 2n+1 polynomials are highlighted in the latter Taylor/interpolation formula: These polynomials can be considered as the shape functions of the numerical technique, these functions taking into account the differential equation via the Taylor formula (5). With the definition (10) (11), the multi-point Taylor series is in the following form: where the polynomials P (x), Q(x) have been computed by the Taylor series of the residual and the Hermite-Birkhoff interpolation. The discretization vector q ∈ R 2n is the unknown at this level and it will be defined in the next part.
2.3. The global problem. In the case of a multi-point Taylor series (n ≥ 2), the two boundary conditions are not sufficient to define the remaining unknowns q ∈ R 2n . Hence 2n − 2 equations are lacking that have to be derived from the differential equation. A secondary cloud of points (y 1 , y 2 , . . . , y 2n−2 ) is introduced where the residual is required to be zero. In summary, the 2n unknowns are obtained by the two boundary conditions and by the differential equation at the 2n−2 points of this secondary cloud: Because of the linearity of the equation (1), the global system (14) is affine with respect to q. The corresponding matrix and right-hand-side vector can be obtained in a straightforward manner. They are given in the Appendix B If necessary, the multi-point Taylor series can be applied in a piecewise manner.
In each interval, a multi-point series is defined as previously. In this case, the global system includes also continuity conditions between neighbor subdomains for the function and its derivative. It can be expressed as: The corresponding matrix and vector are detailed in Appendix B. In this paper, simple pointwise continuity conditions have been considered but other approaches could be used if necessary [2]. In summary, the discrete problem is defined by three families of points, see Figure 1.
• A main cloud of points x 1 , x 2 , . . . , x n where the derivatives of the residual up to order p-2 are vanished, see (5). There are at least two such points in each subdomain.
• A secondary cloud of points where only the residual is vanished.
3. Numerical evaluation. Several numerical tests will be carried out in order to assess the ability of multi-point Taylor series to solve various differential equations. Lopez et al [15] have considered two-point and three-point Taylor series in cases where the one-point series diverges because of the presence of a close singularity in the complex plane. We shall consider here a similar test. Another case will be studied, where the solution is nearly constant in the center of the domain and there are boundary layers. A last example is motivated by the treatment of discontinuities or singularities that can be generated by corners, cracks, pointwise forces. . . To assess the performance of multi-point series in this case, the response to a Dirac distribution will be analyzed numerically. The method of multi-point series offers many variants: single domain or several subdomains, number of expansion points "n", degree "p" of the associated one-point series, location of the clouds of collocation points: these discretization data are the equivalent of the mesh in a finite element method. The numerical tests will permit to discuss the optimal choices and the robustness of the method. The improvement of the discretization is easy within meshless methods but this requires relevant error estimation techniques: in this respect the connection between profile of the residual function (2) and discretization error will be discussed carefully. For simplicity, a symbolic code has been applied in most of the cases: this avoids round-off errors, but can lead to huge computation times for large values of the numbers "n" and "p".  Figure 2 for a basic degree p = 5, 10, 20 and compared with the exact solution in the case p = 5. One sees that the two-point Taylor begins to converge for p = 5 while the two other methods remain wrong. For a large degree p = 20, the three results are very good, but the two-point method provides a very much larger accuracy (10 −9 ) than the one-point series (10 −3 ), the method with two subdomains being in the middle (10 −6 ). Of course the efficiency of the two-point method is certainly due the larger degree of the polynomial (n(p + 1) − 1 instead of p). That is why one can compare the one-point series for p = 20 and the two-point series for p = 10, the degree of the polynomials being about the same (respectively 20 and 21): the two-point series is slightly better in this case (10 −4.7 instead of 10 −3.5 ). This demonstrates the potentialities of the two-point approach, even if the one-point series converges.
In the multi-point approach, one increases the degree of the final polynomial n(p + 1) − 1 by increasing either the number of derivatives "p" or the number "n" of expansion points. For this reason, we plotted the evolution of the error and of the residual R(u(x)) as a function, first of "p" with two points (n = 2), second of the number of expansion points with a constant number of derivatives (p = 4), see    Figure 4. Clearly the convergence is exponential in the two cases. It seems that one converges faster with n than with p: for instance, for n = 5, p = 4, the error is much smaller (10 −9 ) than for n = 2, p = 10 (10 −5 ), the degree of u(x) being about the same (respectively 23 and 21).
One observes also that the maximal residual and the error evolve together in such a way that the residual can be considered as a relevant error indicator. Nevertheless the criterion of the maximal residual has to be handled carefully: indeed with Taylor series, the residual is often very localized and its effect on the solution depends on this localization: a very localized residual has less influence on the solution.
Last the method could be sensitive to the location of the clouds of points because the convergence of a Taylor series depends mainly on the distance between the expansion point and the singularities of the function. Let us consider a case where the approximate solution is not very good (n=2, p=6). The results of Figure 3, 4 have been obtained with a regularly distributed cloud of points x i = ±10/3. The distribution of error and residual for n = 2, p = 6 is represented in Figure 5.The shapes of the two curves are quite similar with a maximal values at about x = ±9. This suggests that the residual is also a good criterion to locate the sources of error. In Figure 6 we have plotted the obtained accuracy as a function of expansion center in the range ±x i ∈ [3.33, 6.5]. An order of magnitude can be gained by passing from a regular distribution x i = ±3.33 to x i = ±5.5. For larger value of x i , the accuracy increases again and the source of error is moved in the center of interval. The best results are found at the transition between maximum error at the end and maximum error in the center.

3.2.
An example where a one-point series is not sufficient. Let us consider the case f (x) = 1, g(x) = 1 This case is referred as Example 2. There is an exact solution u(x) = 1 (x 2 +1) 2 that has two singularities    x = ±i in the complex plane. Because these singularities are close to the domain, a one-point series cannot converge. In this case, our calculations were limited to rather small order (2 ≤ p ≤ 8) because of the limitation due to symbolic computation. The two-point series does not work either and the solution of the three-point series are not completely satisfactory: the error lies between 0, 8% and 24% and the maximal residual is too large, generally very localized near the boundary. Let us consider a four-point series with a small degree and a regularly distributed cloud of expansion   . This is not sufficient to get an accurate solution, the maximal error being 3%. Moreover the maximal residual remains large, but located in narrow regions near the boundary. This is a typical behavior of Taylor series that are accurate close to the expansion point, but may be wrong away. To address this difficulty, the cloud of points was slightly moved to the ends by keeping the degree p = 4. This cloud adaptation has permitted to get acceptable solutions (error= 0.8%, residual= 0.08) with expansion points located in (±0.8, ±1.8). With the latter expansion points, the accuracy decreases exponentially with the degree (Figure 7). Finally one obtains an accuracy less than 10 −3 for p = 6 (degree of the polynomial 27) and less than 10 −4 for p = 8 (degree of the polynomial 35). Thus in view of future 2D/3D applications, the analysis of the residual seems a relevant criterion for an adaptive control, with an adaptivity that can concern the degree p, the number of points n or the location of these points.  3.3. Response to a concentrated force. In the last example (example 3), one discusses the ability of multi-point Taylor series to describe the solution of nonsmooth problems, like the effect of a domain with corners, of cracks or of contact mechanics. More specifically we try to represent the response to a concentrated force, where f (x) = 0, g(x) = δ(x), L = 10, δ(x) being the Dirac distribution. The solution is u(x) = 10−|x| 2 . Such a non-smooth response has to be regularized in order to be represented by an analytic function, as this was done for instance in contact mechanics [1]. The simplest way is to replace the Dirac distribution by a smooth and localized function as δ a (x) = 1 a √ π exp(− x 2 a 2 ), that is known to converge to the Dirac function when the length "a" tends to zero. The solution with the regularized right-hand side g(x) = δ a (x) involves the antiderivative of an error function erf(x) because the latter is the antiderivative of the Gaussian function. Because this solution with g(x) = δ a (x) is not obviously available, the accuracy of multi-point Taylor series is discussed, first with the respect to the solution of the non-smooth problem u(x) = 10−|x| 2 , second by considering the size of the residual. Multi-point Taylor series have been applied in the two cases a = 1, a = 0.1,the last case corresponding to a very localized force. The challenge is to find polynomials approximating these almost discontinuous functions . We were not able to obtain relevant solutions without splitting into subdomains, due to the co-existence of regions of quasi-linearity and regions of very strong variations.   Other data are p = 5, n ∈ [3,5], which leads to degrees of polynomials between 17 and 29. The expansion points x j are regularly distributed in each interval. The functions obtained by these three multi-point series are presented in Figure 8 and 9. One sees that one has been able to represent very localized responses by this procedure. In the case a = 1, one recovers about the same results for n = 3, 4, 5 and a residual always lower than 2.10 −3 (1.310 −5 for n = 4 points by subdomain). The distribution of this residual is plotted in Figure 11 and it appears that the maximal residual is located near the interfaces, i.e. away from the expansion points, as in Example 1. Nevertheless the smallness of the residual and the stability with respect to "n" indicate that the difference with the non-smooth solution is due to the regularization of the r.h.s (a = 0) and not to the approximation by Taylor series. In the second case a = 0.1 corresponding to  the largest localization, the accuracy is acceptable but a bit lower than with a = 1: the maximal residual varies between 7.610 −3 for n = 3 and 2.210 −3 for n = 5. Because of the regularization, the numerical solutions are slightly different from the non-smooth one u(x) = 10−|x| 2 , but they are close: the difference in the middle is less than 6% for a = 1 (u max = 4.718 instead of 5) and less than 1.6% for a = 0.1 (u max = 5.0788).
Likely we have not yet found the best strategies for splitting in subdomains, the best choices for the degree and the number of expansion points, the optimal locations of the discretization points, or the robustness of the method, but this first attempt establishes the ability of the multi-point Taylor method to approximate very localized evolutions.

Conclusion.
A new numerical procedure has been proposed to solve Ordinary Differential Equations and it was evaluated from relevant numerical experiments. The basis is a multi-point Taylor expansion that combines one-point expansion and Hermite-Birkhoff interpolation. It proved effective and reliable. It works well in the simplest cases where a one-point series would be sufficient, but it can converge also exponentially with the degree in cases where any one-point series diverges. It is a meshfree method where the discretization is defined by clouds of points. Hence the implementation is easy and the accuracy can be improved simply by adding new points or by moving these points. In most of the presented examples, there is an exact solution, but it has been established that the error can be estimated and located in other cases, because the residual seems to be a reliable measure of the accuracy and therefore a relevant error estimator. In 3.3, a test involving the regularized Dirac distribution shows that multi-point Taylor expansion is a good candidate to approximate discontinuous functions by high degree polynomials, what could be useful in many applications, for instance in contact mechanics.
This combination of one-point Taylor and high order interpolation could be extended in a multi-dimensional framework. Indeed one-point expansions have been already applied to PDE's [24,22,23] and an interpolation by multivariate polynomials seems well established, see for instance [6,7]. Of course, a full comparison "1-point versus multi-point series" seems very early at this stage, but according to a recent study for 3D PDEs [23], it is clear that the computation time of the Taylor series is a very small part of the full computation time. Always according to [23], a good criterion for such a comparison should be matrix ill-conditioning, but one expects an improved conditioning with multi-point Taylor series.
Appendix A. Hermite-Birkhoff interpolation [18]. The classical interpolation polynomials of Hermite-Birkhoff are recalled for completeness. The data are the degree p and a cloud of n points x 1 , x 2 , . . . , x n . A family of n(p + 1) polynomials permits to interpolate from the values of the function and its derivatives up to order p at this cloud. First one defines a first family of polynomials For instance in the case of two points n = 2, x 1 = − 1 3 , x 2 = 1 3 and of a degree p = 2, the formula (16) yields six polynomials that form a basis for the vector space of polynomials of degree 5: From this first family, one builds the interpolants used in Section 2.2. The last ones (i.e. the first index is maximum) are extracted from (16): and the other follow from a recurrence formula: In the case n = 2, x 1 = −1 3 et x 2 = 1 3 , the six interpolation polynomials are: Appendix B. Details of the global problems. After the local resolution by the method of Taylor series, the remaining equations (14) or (15) are linear with respect to u(x) that is an affine function of the discrete variables q. Hence this leads to a linear problem case K.q = F involving a matrix and a vector. Let us introduce the linear operator associated with the affine operator case : R(u)(x) = − d 2 u(x) dx 2 + f (x)u(x). Thus the discretization of the unknown leads to the discretization of the residuals: R(u) = R(P )(x) + 2n j=1 R(Q j )(x)q j = R(P )(x) + t R(Q)(x).q (22) Let us begin with the case without splitting in subdomains, where the global problem (14) accounts for boundary conditions and collocation inside the domain. With account of (22), this global problem is as follows: . . .
. . . R(P )(y 2(n−1) ) In the case with several subdomains, an approximation in the form (13) is defined in each interval: u k (x) = P k (x) + t Q k (x)q k . The global system (15) includes, first the two boundary conditions and the collocation equations at the secondary clouds of points y 1 , y 2 , . . . , y 2(n−1) , what will be assembled in the matrix K coll and the rhs F coll , second the continuity conditions at the points z k , 1 ≤ k ≤ m − 1, what will be assembled in the matrix K cont and the rhs F cont : . . . P m−1 (z m−1 ) − P m (z m−1 )