# Convergent Power Series ofsech⁡(x) and Solutions to Nonlinear Differential Equations

Authors: U. Al Khawaja; Qasem M. Al-Mdallal
Journal: International Journal of Differential Equations
 (2018)

---

## Abstract

It is known that power series expansion of certain functions such assech⁡(x) diverges beyond a finite radius of convergence. We present here an iterative power series expansion (IPS) to obtain a power series representation of sech⁡(x) that is convergent for all x. The convergent series is a sum of the Taylor series of sech⁡(x) and a complementary series that cancels the divergence of the Taylor series for x≥π/2. The method is general and can be applied to other functions known to have finite radius of convergence, such as 1/(1+x2). A straightforward application of this method is to solve analytically nonlinear differential equations, which we also illustrate here. The method provides also a robust and very efficient numerical algorithm for solving nonlinear differential equations numerically. A detailed comparison with the fourth-order Runge-Kutta method and extensive analysis of the behavior of the error and CPU time are performed.


---

## Body

## 1. Introduction

It is well-known that the Taylor series of some functions diverge beyond a finite radius of convergence [1]. For instance, by way of example not exhaustive enumeration, the Taylor series of sech⁡(x) and 1/(1+x2) diverge for x≥π/2 and x≥1, respectively. Increasing the number of terms in the power series does not increase the radius of convergence; it only makes the divergence sharper. The radius of convergence can be increased only slightly via some functional transforms [2]. Among the many different methods of solving nonlinear differential equations [3–9], the power series is the most straightforward and efficient [10]. It has been used as a powerful numerical scheme for many problems [11–19] including chaotic systems [20–23]. Many numerical algorithms and codes have been developed based on this method [10–12, 20–24]. However, the above-mentioned finiteness of radius of convergence is a serious problem that hinders the use of this method to wide class of differential equations, in particular the nonlinear ones. For instance, the nonlinear Schrödinger equation (NLSE) with cubic nonlinearity has the sech⁡(x) as a solution. Using the power series method to solve this equation produces the power series of a sech⁡(x), which is valid only for x<π/2.A review of the literature reveals that the power series expansion was exploited by several researchers [10–12, 20–24] to develop powerful numerical methods for solving nonlinear differential equations. Therefore, this paper is motivated by a desire to extend these attempts to a develop a numerical scheme with systematic control on the accuracy and error. Specifically, two main advances are presented in this paper: (1) a method of constructing a convergent power series representation of a given function with an arbitrarily large radius of convergence and (2) a method of obtaining analytic power series solution of a given nonlinear differential equation that is free from the finite radius of convergence. Through this paper, we show robustness and efficiency of the method via a number of examples including the chaotic Lorenz system [25] and the NLSE. Therefore, solving the problem of finite radius of convergence will open the door wide for applying the power series method to much larger class of differential equations, particularly the nonlinear ones.It is worth mentioning that the literature includes several semianalytical methods for solving nonlinear differential equations; such as homotopy analysis method (HAM), homotopy perturbation method (HPM), and Adomian decomposition method (ADM); for more details see [26–29] and the references therein. Essentially, these methods generate iteratively a series solution for the nonlinear systems where we have to solve a linear differential equation at each iteration. Although these methods prove to be effective in solvingmost of nonlinear differential equations and in obtaining a convergent series solution, they have few disadvantages such as the large number of terms in the solution as the number of iterations increases. One of the most important advantages of the present technique is the simplicity in transforming the nonlinear differential equation into a set of simple algebraic difference equations which can be easily solved.The paper is thus divided into two, seemingly separated, but actually connected main parts. In the first (Section2), we show, for a given function, how a convergent power series is constructed out of the nonconverging one. In the second part (Section 3.1), we essentially use this idea to solve nonlinear differential equations. In Section 3.2, we investigate the robustness and efficiency of the method by studying the behavior of its error and CPU time versus the parameters of the method. We summarise our results in Section 4.


## 2. Iterative Power Series Method

This section describes how to obtain a convergent power series for a given function that is otherwise not converging for allx. In brief, the method is described as follows. We expand the function f(x) in a power series as usual, say around x=0. Then we reexpress the coefficients, f(n)(x), in terms of f(x). This establishes a recursion relation between the higher-order coefficients, f(n)(0), and the lowest order ones, f(0)(0) and f(1)(0), and thus the power series is written in terms of only these two coefficients. Then the series and its derivative are calculated at x=Δ, where Δ is much less than the radius of convergence of the power series. A new power series expansion of f(x) is then performed at x=Δ. Similarly, the higher-order coefficients are reexpressed in terms of the lowest order coefficients f(0)(Δ) and f(1)(Δ). The value of the previous series and its derivative calculated at x=Δ are then given to f(0)(Δ) and f(1)(Δ), respectively. Then a new expansion around 2Δ is performed with the lowest order coefficients being taken from the previous series, and so on. This iterative process is repeated N times. The final series will correspond to a convergent series at x=NΔ.Here is a detailed description of the method. The functionf(x) is expanded in a Taylor series, T0(x), around x=0. The infinite Taylor series is an exact representation of f(x) for x<R where R is the radius of convergence. For x≥R the series diverges. We assume that x is divided into N small intervals Δ=x/N such that Δ<R. Expanding f(x) around the beginning of each interval we obtain N convergent Taylor series representing f(x) in each interval(1)Tjy=∑n=0∞1n!fnjΔy-jΔn,jΔ≤y<j+1Δ,j=0,1,2,…,N,where Tj(y) denotes the Taylor series expansion of f(y) around y=jΔ and f(n)(jΔ) is the nth derivative of f(y) calculated at y=jΔ. It is noted that we use y∈[(j-1)Δ,jΔ] as the independent variable for the nth Taylor series expansion to distinguish it from x=NΔ. However, these series can not be combined in a single series since their ranges of applicability are different and do not overlap. To obtain a single convergent power series out of the set of series Tj, we put forward two new ideas, which constitute the basis of our method; namely:(
1
) Reexpress f(n)(y) in terms of f(y) as f(n)(y)=Fn[f(y)], where the functional Fn[f(y)] is determined by direct differentiation of f(y) for n times and then reexpressing the result in terms of f(y) only. We conjecture that this is possible for a wide class of functions if not all. At least for the two specific functions considered here, this turned out to be possible. Equation (1) then takes the form(2)Tjy=∑n=0∞ana0jy-jΔn,where we have renamed f(jΔ) by a0j and Fn[f(jΔ)]/n! by an(a0j) for a reason to be obvious in the next section. Thus, the coefficients an for all n are determined only by a0j.(
2
) Calculate a0j from Tj-1 at jΔ(3)a0j=Tj-1jΔ=∑n=0∞ana0j-1Δn,j=1,2,…,N,which amounts to assigning the value of the Taylor series at the end of an interval to a0j of the consecutive one. Equation (3) captures the essence of the recursive feature of our method; a0N is calculated recursively from a00 by repeated action of the right-hand-side on a00. While Tj represents the function f within an interval of width Δ, the sequence a0j corresponds to the values of the function at the end of the intervals. In the limit N→∞, or equivalently Δ→0, the discrete set of a0j values and jΔ render to the continuous function f(x) and its independent variable x, respectively. Formally, the convergent power series expansion of f(x) around x=0 will thus be given by(4)fx=limN→∞⁡SN,where SN denotes the Nth iteration of(5)Sf0=∑n=0∞anf0xNn.As an illustrative example, we apply the method tof(x)=sech⁡(x). The infinite Taylor series expansion of this function diverges sharply to infinity at x=π/2. The first step is to determine the coefficients anj, which are the coefficients of the Tj series(6)Tjy=sech⁡yy=jΔ+-sech⁡ytanh⁡yy=jΔy-jΔ+12!-sech3⁡y+sech⁡ytanh2⁡yy=jΔy-jΔ2+13!5sech2⁡y-tanh2⁡ysech⁡ytanh⁡yy=jΔy-jΔ3+14!5sech5⁡y-18sech3⁡ytanh2⁡y+sech⁡ytanh4⁡yy=jΔy-jΔ4+⋯.The next step is to reexpress the higher-order coefficients, anj, in terms of the zeroth-order coefficient a0j=sech⁡(y)y=jΔ. The property sech2⁡y+tanh2⁡y=1 is used to that end. It is noticed, however, that while it is possible to express the even-n coefficients in terms of a0j only, the odd-n coefficients can only be expressed terms of both a0j and a1j=-sech⁡ytanh⁡(y)y=jΔ. In the context of solving differential equations using the power series method, this reflects the fact that the solution is expressed in terms of two independent parameters (initial conditions). The sech function is indeed a solution of a second-order differential equation, which is solved using this method in the next section. Equation (6) then takes the form(7)Tjy=a0j+a1jy-jΔ+12!1-2a0j2a0jy-jΔ2+13!1-6a0j2a1jy-jΔ3+14!1-8a0j2+8a0j4-12a1j2a0jy-jΔ4+⋯.Calculating Tj(y) series at the end of its interval of applicability, y=(j+1)Δ, we get(8)Tjj+1Δ=a0j+a1jΔ+a2ja0jΔ2+a3ja0j,a1jΔ3+a4ja0jΔ4+⋯,where the “recursion” coefficients are given by(9)a2ja0j=12!1-2a0j2a0j,a3ja0j,a1j=13!1-6a0j2a1j,a4ja0j=14!1-8a0j2+8a0j4-12a1j2a0j.Finally, we assign Tj((j+1)Δ) to a0j+1(10)a0j+1=a0j+a1jΔ+a2ja0jΔ2+a3ja0j,a1jΔ3+a4ja0jΔ4+⋯.The second independent coefficient a1j+1 is determined by the derivative of Tj(y) calculated at y=(j+1)Δ(11)a1j+1=a1j+2a2ja0jΔ+3a3ja0j,a1jΔ2+4a4ja0jΔ3+⋯.While, in the limit N→∞, a0j corresponds to the function f(x), the sequence a1j corresponds to f′(x). Therefore, the power series expansion of sech⁡(x) and its first derivative are given by(12)fxf′x=limN→∞⁡a0+a1xN+a2a0xN2+a3a0,a1xN3+a4a0xN4+⋯a1+2a2a0xN+3a3a0,a1xN2+4a4a0xN3+⋯N,where the superscript of the matrix on the right-hand-side, N, denotes the Nth iteration of the matrix. The superscript j has been removed since the functional form of the recursion coefficients does not depend on j. The procedure of calculating the power series recursively is described as follows. First, a0=sech⁡(0)=1 and a1=sech′⁡(0)=0 are substituted in the right-hand-side of the last equation. Then the result of the upper element is taken as the updated value of a0, and, similarly, the lower element updates a1. The two updated values are then resubstituted back in the right-hand-side. The process is repeated N times. To obtain an explicit form of the series we truncate the Taylor series at nmax=4 and use N=4 iterations. The resulting expansion takes the form(13)sech⁡x=1-12x2+524x4-0.0806681x6+0.0302048x8+⋯+1.4434798×10-461x624.It is noted that the higher-order coefficients, which correspond to ratios of large integers, are represented in real numbers for convenience. Already with such a small number of iterations, N=4, the number of terms equals 313. By inspection, we find that the number of terms equals ((nmax+1)N+1)/2. Here, nmax is even due to the fact that sech⁡(x) is an even function.It is also noted that the series (13) is composed of the Taylor expansion of sech⁡(x) around zero, represented by the first three terms, and a series of higher-order terms generated from the nonlinearity in the recursion relations of an. In fact, we prove in the next section that this property holds for any nmax, N, and function f(x), provided that the Taylor series of the later exists. Therefore, the power series expansion of sech⁡(x), given by (12), can be put in the suggestive form(14)sech⁡x=T+limN→∞⁡CN,where T is the infinite Taylor series of f(x) about x=0 and C(N) is a complementary series. It turns out that the complementary series increases the radius of convergence of T for x≥π/2. For finite N, the effect of C(N) is to shift the radius of convergence, R, to a larger value such that R→∞ for N→∞. In Figure 1 we plot the convergent power series obtained by the present method as given by (12) using nmax=4 and N=100. The curve is indistinguishable from the plot of sech⁡(x). Both the Taylor series expansion, T, and the complementary series, C, diverge sharply at x=π/2. Since C is essentially zero for x<π/2, it will not affect the sum T+C. However, its major role is to cancel the divergency for x≥π/2. In the limit N→∞, T will be an exact representative of sech⁡(x) for x<π/2 and C will equal zero in the same interval. For x≥π/2, the divergences in T and C cancel each other with a remainder that is an exact representative of sech⁡(x). In this manner, T+C will represent sech⁡(x) for all x.Figure 1
The solid curve corresponds to the convergent power series obtained by the present method as given by (12) using nmax=4 and N=100. It is indistinguishable from the solid curve of f(x)=sech⁡(x). Dashed curve corresponds to the Taylor series expansion, T, and the dotted curve corresponds to the complementary series, C.For finite values ofnmax and N, the series T+C is an approximate representative of sech⁡(x). Truncating the Taylor series at nmax introduces an error of order Δnmax+1. This error will be magnified N times due the recursive substitutions. The total error is then estimated by(15)error=xNnmax+1N.For the parameters used in Figure 1, this error is of order 10-6 at x=5. This can be reduced to extremely small values such as 10-131 with nmax=100. However, the number of terms in the series T+C will be of order 10200 which is extremely large and hinders any analytical manipulations.As another example, we considerf(x)=1/(1+x2) with Taylor series diverging at x=1. Much of the formulation we followed for the previous case holds here and the specifics of the function alter only the recursion relations, (9):(16)a2ja0j=3-4a0ja0j2,(17)a3ja0j,a1j=21-2a0ja0ja1j,(18)a4ja0j=16-64a0j+85a0j2-52a0j3+16a0j4a0j.The convergent power series is obtained by using these recursion relations in (12). Plots similar to those of Figure 1 are obtained.We present now a proof that the convergent power series produced by the recursive procedure always regenerates the Taylor series in addition to a complementary one.Proposition 1.
If we expandf(x) in a Taylor series, T, around x=0 truncated at nmax and use the recursive procedure, as described above, the resulting convergent power series always takes the form T+C where C is a power series of orders larger than nmax. This is true for any number of iterations, N, maximum power of the Taylor series, nmax, and for all functions that satisfy the general differential equation(19)f′′x=Ffx,where F[·] is an analytic real functional that does not contain derivatives.Proof.
It is trivial to prove this for a specific case, such assech⁡(x). For the general case, we prove this only for nmax=4 and N=2. The Taylor series expansion of sech⁡(x) around x=0 is(20)sech⁡x=1-12x2+524x4+⋯.In our recursive procedure, this is put in the equivalent form(21)sech⁡Δsech′⁡Δ≃a0+a1Δ+1-2a022Δ2+13!a11-6a02Δ3+14!a01-8a02+12a04-12a12Δ41-2a02Δ+12a11-6a02Δ2+16a01-8a02+12a04-12a12Δ3N,where the approximation stems from using finite N and nmax, and Δ=x/N. For N=1, a0=1, and a1=0, (20) is regenerated. However, in our recursive procedure a0 and a1 are kept as variables since they will be substituted for at each recursive step. Only at the last step are their numerical values inserted. For N=2, we resubstitute in the last equation for a0 and a1 by their updated expressions, as follows:(22)a0a1⟶a0+a1Δ+1-2a022Δ2+13!a11-6a02Δ3+14!a01-8a02+12a04-12a12Δ41-2a02Δ+12a11-6a02Δ2+16a01-8a02+12a04-12a12Δ3.Substituting the updated expressions for a0 and a1 in (21), we get(23)sech⁡x≃a0+a1xN+2a01-2a02xN2+43a11-6a02xN3+23a01-8a02+12a04-12a12xN4.Clearly for N=2, the last equation gives the Taylor expansion, that is, (21) with N=1. The complimentary series, C, is absent here since we have terminated the expansions at n=nmax=4. For N=3, another step of resubstituting updated expressions is needed, and so on.
Now, we present the proof for the more general case, namely, whenf(x) is unspecified but is a solution to (19). We start with the following Taylor series expansion of f(x) and its derivative(24)fΔf′Δ=a0+a1Δ+a2a0,a1Δ2+a3a0,a1Δ3+a4a0,a1Δ4a1+2a2a0,a1Δ+3a3a0,a1Δ2+4a4a0,a1Δ3.Substituting on the right-hand-side for a0 and a1 by f(Δ) and f′(Δ), respectively, we get(25)fΔ=a0+a1Δ+4a2Δ2+5a3+2a2∂a2∂a1+a1∂a2∂a0Δ3+6a4+3a3∂a2∂a1+2a2∂a3∂a1+2a22∂2a2∂a12+a2∂a2∂a0+a1∂a3∂a0+2a1a2∂2a2∂a1∂a0+12a12∂2a2∂a02Δ4.The partial derivatives can not be calculated unless the functional forms of the recursion coefficients are known. One possibility is to specify the function being expanded, f(x), as we did at the start of this proof. The other possibility is to exploit the differential equation that f(x) is a solution for, namely, (19). Substituting the Taylor expansion of f(Δ), from (24) in (19) and expanding up to the fourth power in Δ, we obtain the following relations:(26)a2=-F2,a3=-a1F′6,a4=-a2F′12-a12F′′24,which lead to(27)∂a2∂a0=-F′2,∂a2∂a1=-ΔF′2,∂2a2∂a02=-F′′2,∂2a2∂a12=-Δ2F′′2,∂2a2∂a0∂a1=-ΔF′′2,∂a3∂a1=-F′6.Using these equations in (28), we obtain(28)fΔ=a0+a1Δ+4a2Δ2+8a3a0,a1Δ3+16a4a0,a1Δ4.For N=2, the last equation regenerates the Taylor series of f(x), namely, (24) with N=1, and this completes the proof.


## 3. Application to Nonlinear Differential Equations

The method described in the previous section can be used as a powerful solver and integrator of nonlinear differential equations both analytically and numerically. In Section3.1, we apply the method on a number of well-known problems. In Section 3.2, we show the power of the method in terms of detailed analysis of the error and CPU time.


### 3.1. Examples

For the sake of demonstration, we consider the following nonlinear differential equation:(29)12fx-12f′′x-f3x=0.The reason for selecting this equation is that f(x)=sech⁡(x) is one of its exact solutions. Substituting the power series expansion f(x)=∑n=0∞anxn, we obtain, as usual, the recursion relations(30)a2a0=12!1-2a02a0,a3a0,a1=13!1-6a02a1,a4a0=14!1-8a02+12a04-12a12a0,where a0 and a1 turn out to be independent parameters which in the present case correspond to the initial conditions on the solution and its first derivative. It is not surprising that these recursion relations are identical with those we found for the sech⁡(x) in the previous section, (9). Therefore, substituting the above recursion relations in f(x)=∑n=0∞anxn we obtain the Taylor series expansion, T, of f(x)=sech⁡(x). Removing the divergency in T follows exactly the same steps as in the previous section, and thus an exact solution in terms of a convergent power series is obtained, as also plotted in Figure 1.Forf(x)=1/(1+x2), the relevant differential equation is(31)f′′x+f3x-6f2x=0.Substituting the power series expansion in this equation, the recursion relations will be given by (16)–(18). Similarly, the convergent series solution will be obtained, as in the previous section.


### 3.2. Numerical Method

As a numerical method, the power series is very powerful and efficient [10]. The power series method with Nmax=4, denoted by IPS4, is used to solve the NLSE, (29) and the error is calculated as the difference between the numerical solution and the exact solution, namely, sech⁡(x). The equation is then resolved using the fourth-order Runge-Kutta (RK4) method. In Figure 2, we plot the error of both methods which turn out to be of the same order. Using the iterative power series method with Nmax=12, (IPS12), the error drops to infinitesimally low values. Neither the CPU time nor the memory requirements for IPS12 are much larger than those for IPS4; it is straight forward upgrade to higher orders which leads to ultrahigh efficiency. This is verified by the set of tables, Tables 1–4, where we compute sech⁡(1) using both the RK4 and the iterative power series method and show the corresponding CPU time. For the same N, Tables 1 and 3 show that both RK4 and IPS4 produce the first 16 digits of the exact value (underlined numbers in the last raw) and consume almost the same CPU time. Table 3 shows that, for the same N, IPS12, reproduces the first 49 digits of the exact value. The CPU time needed for such ultrahigh accuracy is just about 3 times that of the RK4 and IPS4. Of course the accuracy can be arbitrarily increased by increasing N or more efficiently Nmax. For IPS12 to produce only the first 16 digits, as in RK4 and IPS4, only very small number of iterations is needed, as shown in Table 4. The CPU time in this case is about 100 times less than that of RK4 and IPS4, highlighting the high efficiency of the power series method.Table 1
RK4sech⁡(x) solution of the NLSE, (29), computed at x=1.
RK4
N
sech
⁡
(
1
)
CPU time
1000
0.6480542736639463753683598987682775440961395993015
0.087257
2000
0.6480542736638892102461002212645001242748400651983
0.166534
3000
0.6480542736638861522791729310382098086291071707752
0.229016
4000
0.6480542736638856377319630763886301512099763161750
0.314480
5000
0.6480542736638854971232600221160524143929321411426
0.427610
Exact
0.6480542736638853995749773532261503231084893120719Table 2
IPS4sech⁡(x) solution of the NLSE, (29), computed at x=1.
IPS4
N
sech
⁡
(
1
)
CPU time
1000
0.6480542736639323955233350367786700400715255548373
0.056793
2000
0.6480542736638883277492809223389308596155602276884
0.137941
3000
0.6480542736638859773823119964459740185046024606302
0.215395
4000
0.6480542736638855823023023050565759829475356280281
0.304501
5000
0.6480542736638854743968579026604641270130595226186
0.372284
Exact
0.6480542736638853995749773532261503231084893120719Table 3
IPS12sech⁡(x) solution of the NLSE, (29), computed at x=1.
IPS12
N
sech
⁡
(
1
)
CPU time
1000
0.6480542736638853995749773532261503231079594354079
0.280617
2000
0.6480542736638853995749773532261503231084891816361
0.595639
3000
0.6480542736638853995749773532261503231084893110634
0.891370
4000
0.6480542736638853995749773532261503231084893120400
1.080357
5000
0.6480542736638853995749773532261503231084893120697
1.366386
Exact
0.6480542736638853995749773532261503231084893120719Table 4
IPS12sech⁡(x) solution of the NLSE, (29), computed at x=1, but with much less number of iterations than in Table 3.
IPS12
N
sech
⁡
(
1
)
CPU time
3
0.6480542794079665629469114154348980055814088430953
0.000782
6
0.6480542736643770346283969779587807646256058100135
0.001429
9
0.6480542736638872007452856567074697922787489590238
0.002123
12
0.6480542736638854259793573015747954030451932565027
0.002768
15
0.6480542736638854001495023257702479909741369016589
0.003594
Exact
0.6480542736638853995749773532261503231084893120719Figure 2
Error defined by the difference between the numerical and exact soliton solution of the NLSE, (29). (a) Dashed line corresponds to RK4 and dotted line corresponds to IPS4. (b) IPS12. Parameters used: N=1000, number of digits Nd=50, f(0)=1, and f′(0)=0.
(a)
(b)A more challenging test on the power series method is the chaotic Lorenz system [25] given by (32)z˙1=-σz1+σz2,z˙2=-z1z3+rz1-z2,z˙3=z1z2-bz3,where we take the usual values σ=10, b=-8/3, and R=28 with initial conditions z1(0)=z2(0)=1 and z3(0)=20. It is straight forward to generalise the method to three differential equations; therefore we do not show the details of the calculation. In Figure 3, the results of solving the Lorenz system using RK4 and IPS12 are shown. For the same parameters, namely discretization, RK4 reaches stability at about x=30; that is, the curve for x<30 is unchanged by increasing N, but for x>30, the curve keeps changing by increasing N. In comparison, IPS12 reaches stability at about x=50. In chaotic systems, it is quite challenging to go that deep in the chaotic region. Hence, the need for such high accuracy methods.Figure 3
Numerical solution of the chaotic Lorenz system, (32), using the RK4 and IPS methods. Parameters: for IPS (solid black line), Nmax=16, N=5000. For RK4 (dashed red line), N=2×105. For both curves we used Nd=1000, z1=1=z2=1, z3=20, σ=10, b=-8/3, and R=28.Achieving higher accuracy requires larger CPU time usage. Therefore, it is important to investigate how the CPU time, denoted here byT, depends on the main parameters of the method, namely N and Nmax. A typical plot is shown in Figure 4, where we plot on a log-log scale the CPU time versus the error. The linear relationship indicates T∝errorp, where p is the slope of the line joining the points in Figure 4. The error can be calculated in two ways: (i) the difference between the numerical solution and (ii) theoretical estimate, (15). Both ways are shown in the figure and they have the same slope. However, as expected, error defined by (15), which is actually an upper limit on the error, is always larger than the first one. To find how the CPU time depends explicitly on N and Nmax, we argue that the dependence should be of the form T∝NNmax3. This is justified by the fact that CPU time should be linearly proportional to the number of terms computed. The number of terms computed increases linearly with the number of iterations N. The number of terms in the power series is linearly proportional to Nmax. When substituted in the NLSE with cubic nonlinearity, the resulting number of terms, and thus T, will be proportional to Nmax3. In Figure 5, it is shown indeed that the ratio T/NNmax3 saturates asymptotically to a constant for large N and Nmax since the scaling behaviors mentioned here apply for large N and Nmax. The proportionality constant, c, is very small and corresponds to the CPU time of calculating one term. It is dependent on the machine, the programming optimization [10], and the number of digits used, Nd. In terms of the number of digits, the CPU time increases, as shown in Figure 6, where it is noticed that CPU time is almost constant for number of digits Nd<500.Figure 4
CPU time versus error on a log-log scale for the iterative power series solution of the NLSE, (29). Black filled circles correspond to error calculated from the difference between the power series solution and the exact solution. Red empty circles correspond to the error estimated by the formula (15). Parameters used: Nmax=50, nd=500, and N ranges between 2 and 90.Figure 5
CPU time versusNmax (a) and N (b) for the NLSE, (29). Parameters used: for (a), N=400. For (b), Nmax=70. For both plots, Nd=1000, c=6.2×10-7, f(0)=1, f′(0)=0, and Δ=1/N.
(a)
(b)Figure 6
CPU time versus the number of digits for the power series solution of NLSE, (29). Parameters used: Nmax=4, N=400, f(0)=1, f′(0)=0, and Δ=1/N. Inset shows a zoom on linear scale. Red dashed line is an asymptote of slope equal 1.


## 3.1. Examples

For the sake of demonstration, we consider the following nonlinear differential equation:(29)12fx-12f′′x-f3x=0.The reason for selecting this equation is that f(x)=sech⁡(x) is one of its exact solutions. Substituting the power series expansion f(x)=∑n=0∞anxn, we obtain, as usual, the recursion relations(30)a2a0=12!1-2a02a0,a3a0,a1=13!1-6a02a1,a4a0=14!1-8a02+12a04-12a12a0,where a0 and a1 turn out to be independent parameters which in the present case correspond to the initial conditions on the solution and its first derivative. It is not surprising that these recursion relations are identical with those we found for the sech⁡(x) in the previous section, (9). Therefore, substituting the above recursion relations in f(x)=∑n=0∞anxn we obtain the Taylor series expansion, T, of f(x)=sech⁡(x). Removing the divergency in T follows exactly the same steps as in the previous section, and thus an exact solution in terms of a convergent power series is obtained, as also plotted in Figure 1.Forf(x)=1/(1+x2), the relevant differential equation is(31)f′′x+f3x-6f2x=0.Substituting the power series expansion in this equation, the recursion relations will be given by (16)–(18). Similarly, the convergent series solution will be obtained, as in the previous section.


## 3.2. Numerical Method

As a numerical method, the power series is very powerful and efficient [10]. The power series method with Nmax=4, denoted by IPS4, is used to solve the NLSE, (29) and the error is calculated as the difference between the numerical solution and the exact solution, namely, sech⁡(x). The equation is then resolved using the fourth-order Runge-Kutta (RK4) method. In Figure 2, we plot the error of both methods which turn out to be of the same order. Using the iterative power series method with Nmax=12, (IPS12), the error drops to infinitesimally low values. Neither the CPU time nor the memory requirements for IPS12 are much larger than those for IPS4; it is straight forward upgrade to higher orders which leads to ultrahigh efficiency. This is verified by the set of tables, Tables 1–4, where we compute sech⁡(1) using both the RK4 and the iterative power series method and show the corresponding CPU time. For the same N, Tables 1 and 3 show that both RK4 and IPS4 produce the first 16 digits of the exact value (underlined numbers in the last raw) and consume almost the same CPU time. Table 3 shows that, for the same N, IPS12, reproduces the first 49 digits of the exact value. The CPU time needed for such ultrahigh accuracy is just about 3 times that of the RK4 and IPS4. Of course the accuracy can be arbitrarily increased by increasing N or more efficiently Nmax. For IPS12 to produce only the first 16 digits, as in RK4 and IPS4, only very small number of iterations is needed, as shown in Table 4. The CPU time in this case is about 100 times less than that of RK4 and IPS4, highlighting the high efficiency of the power series method.Table 1
RK4sech⁡(x) solution of the NLSE, (29), computed at x=1.
RK4
N
sech
⁡
(
1
)
CPU time
1000
0.6480542736639463753683598987682775440961395993015
0.087257
2000
0.6480542736638892102461002212645001242748400651983
0.166534
3000
0.6480542736638861522791729310382098086291071707752
0.229016
4000
0.6480542736638856377319630763886301512099763161750
0.314480
5000
0.6480542736638854971232600221160524143929321411426
0.427610
Exact
0.6480542736638853995749773532261503231084893120719Table 2
IPS4sech⁡(x) solution of the NLSE, (29), computed at x=1.
IPS4
N
sech
⁡
(
1
)
CPU time
1000
0.6480542736639323955233350367786700400715255548373
0.056793
2000
0.6480542736638883277492809223389308596155602276884
0.137941
3000
0.6480542736638859773823119964459740185046024606302
0.215395
4000
0.6480542736638855823023023050565759829475356280281
0.304501
5000
0.6480542736638854743968579026604641270130595226186
0.372284
Exact
0.6480542736638853995749773532261503231084893120719Table 3
IPS12sech⁡(x) solution of the NLSE, (29), computed at x=1.
IPS12
N
sech
⁡
(
1
)
CPU time
1000
0.6480542736638853995749773532261503231079594354079
0.280617
2000
0.6480542736638853995749773532261503231084891816361
0.595639
3000
0.6480542736638853995749773532261503231084893110634
0.891370
4000
0.6480542736638853995749773532261503231084893120400
1.080357
5000
0.6480542736638853995749773532261503231084893120697
1.366386
Exact
0.6480542736638853995749773532261503231084893120719Table 4
IPS12sech⁡(x) solution of the NLSE, (29), computed at x=1, but with much less number of iterations than in Table 3.
IPS12
N
sech
⁡
(
1
)
CPU time
3
0.6480542794079665629469114154348980055814088430953
0.000782
6
0.6480542736643770346283969779587807646256058100135
0.001429
9
0.6480542736638872007452856567074697922787489590238
0.002123
12
0.6480542736638854259793573015747954030451932565027
0.002768
15
0.6480542736638854001495023257702479909741369016589
0.003594
Exact
0.6480542736638853995749773532261503231084893120719Figure 2
Error defined by the difference between the numerical and exact soliton solution of the NLSE, (29). (a) Dashed line corresponds to RK4 and dotted line corresponds to IPS4. (b) IPS12. Parameters used: N=1000, number of digits Nd=50, f(0)=1, and f′(0)=0.
(a)
(b)A more challenging test on the power series method is the chaotic Lorenz system [25] given by (32)z˙1=-σz1+σz2,z˙2=-z1z3+rz1-z2,z˙3=z1z2-bz3,where we take the usual values σ=10, b=-8/3, and R=28 with initial conditions z1(0)=z2(0)=1 and z3(0)=20. It is straight forward to generalise the method to three differential equations; therefore we do not show the details of the calculation. In Figure 3, the results of solving the Lorenz system using RK4 and IPS12 are shown. For the same parameters, namely discretization, RK4 reaches stability at about x=30; that is, the curve for x<30 is unchanged by increasing N, but for x>30, the curve keeps changing by increasing N. In comparison, IPS12 reaches stability at about x=50. In chaotic systems, it is quite challenging to go that deep in the chaotic region. Hence, the need for such high accuracy methods.Figure 3
Numerical solution of the chaotic Lorenz system, (32), using the RK4 and IPS methods. Parameters: for IPS (solid black line), Nmax=16, N=5000. For RK4 (dashed red line), N=2×105. For both curves we used Nd=1000, z1=1=z2=1, z3=20, σ=10, b=-8/3, and R=28.Achieving higher accuracy requires larger CPU time usage. Therefore, it is important to investigate how the CPU time, denoted here byT, depends on the main parameters of the method, namely N and Nmax. A typical plot is shown in Figure 4, where we plot on a log-log scale the CPU time versus the error. The linear relationship indicates T∝errorp, where p is the slope of the line joining the points in Figure 4. The error can be calculated in two ways: (i) the difference between the numerical solution and (ii) theoretical estimate, (15). Both ways are shown in the figure and they have the same slope. However, as expected, error defined by (15), which is actually an upper limit on the error, is always larger than the first one. To find how the CPU time depends explicitly on N and Nmax, we argue that the dependence should be of the form T∝NNmax3. This is justified by the fact that CPU time should be linearly proportional to the number of terms computed. The number of terms computed increases linearly with the number of iterations N. The number of terms in the power series is linearly proportional to Nmax. When substituted in the NLSE with cubic nonlinearity, the resulting number of terms, and thus T, will be proportional to Nmax3. In Figure 5, it is shown indeed that the ratio T/NNmax3 saturates asymptotically to a constant for large N and Nmax since the scaling behaviors mentioned here apply for large N and Nmax. The proportionality constant, c, is very small and corresponds to the CPU time of calculating one term. It is dependent on the machine, the programming optimization [10], and the number of digits used, Nd. In terms of the number of digits, the CPU time increases, as shown in Figure 6, where it is noticed that CPU time is almost constant for number of digits Nd<500.Figure 4
CPU time versus error on a log-log scale for the iterative power series solution of the NLSE, (29). Black filled circles correspond to error calculated from the difference between the power series solution and the exact solution. Red empty circles correspond to the error estimated by the formula (15). Parameters used: Nmax=50, nd=500, and N ranges between 2 and 90.Figure 5
CPU time versusNmax (a) and N (b) for the NLSE, (29). Parameters used: for (a), N=400. For (b), Nmax=70. For both plots, Nd=1000, c=6.2×10-7, f(0)=1, f′(0)=0, and Δ=1/N.
(a)
(b)Figure 6
CPU time versus the number of digits for the power series solution of NLSE, (29). Parameters used: Nmax=4, N=400, f(0)=1, f′(0)=0, and Δ=1/N. Inset shows a zoom on linear scale. Red dashed line is an asymptote of slope equal 1.


## 4. Conclusions

We have presented an iterative power series method that solves the problem of finite radius of convergence. We have proved that the iterative power series is always composed of a sum of the typical power series of the function and a complementary series that cancels the divergency. The method is divided into two schemes where in the first we find a convergent power series for a given function and in the second we solve a given nonlinear differential equation. The result of the iterative power series expansion ofsech⁡(x) is remarkably convergent for arbitrary radius of convergence and accuracy, as shown by Figures 1 and 2 and Tables 1–4. Extremely high accuracy can be obtained by using higher-order iterative power series via increasing Nmax with relatively low CPU time usage. Robustness and efficiency of the method have been shown by solving the chaotic Lorenz system and the NLSE. Extensive analysis of the error and CPU time characterising the method is performed. Although we have focused on the localised sech⁡(x) solution of the NLSE, all other solitary wave solutions (conoidal waves) can be obtained using the present method, just by choosing the appropriate initial conditions.The method can be generalised to partial and fractional differential equations making its domain of applicability even wider.



---

*Source: 6043936-2018-02-13.xml*