Weak periodic solutions of $x\ddot{x} + 1 = 0$ and the harmonic balance method

We prove that the differential equation $x\ddot{x} + 1 = 0$ has continuous weak periodic solutions and compute their periods. Then, we use the Harmonic Balance Method until order six to approach these periods and to illustrate how the sharpness of the method increases with the order. Our computations rely on the Gr\"obner basis method.


Introduction and main results
The nonlinear differential equation appears in the modeling of certain phenomena in plasma physics [1]. In [6], Mickens calculates the period of its periodic orbits and also uses the N-th order Harmonic Balance Method (HBM), for N = 1, 2, to obtain approximations of these periodic solutions and of their corresponding periods. Strictly speaking, it can be easily seen that neither equation (1), nor its associated system ẋ = y, which is singular at x = 0, have periodic solutions. Our first result gives two different interpretations of Mickens' computation of the period. The first one in terms of weak (or generalized) solutions. In this work a weak solution will be a function satisfying the differential equation (1) on an open and dense set, but being of class C 0 at some isolated points. The second one, as the limit, when k tends to zero, of the period of actual periodic solutions of the extended planar differential system which, for k = 0, has a global center at the origin. Recall that the N-th order HBM consists in approximating the solutions of differential equations by truncated Fourier series with N harmonics and an unknown frequency; see for instance [7,8] or Section 3 for a short overview of the method. In [8, p. 180] the author asks for techniques for dealing analytically with the N-th order HBM, for N ≥ 3. In [4] it is shown how resultants can be used when N = 3. Here we utilize a more powerful tool, the computation of Gröbner basis ([2, Ch. 5]), for going further in the obtention of approximations of the function T (A) introduced in Theorem 1.1.
Notice that equation (1) is equivalent to the family of differential equations for any m ∈ N ∪ {0}. Hence it is natural to approach the period function, by the periods of the trigonometric polynomials obtained applying the N-th order HBM to (4). Next theorem gives our results for N ≤ 6. Here [a] denotes the integer part of a.
(ii) For m = 0, (iv) For m = 2, Moreover, the approximate values appearing above are roots of given polynomials with integer coefficients. Whereby the Sturm sequences approach can be used to get them with any desirable precision.
Notice that the values T 1 (A; m), for m ∈ {0, 1, 2} given in items (ii), (iii) and (iv), respectively, are already computed in item (i). We only explicite them to clarify the reading.
Observe that the comparison of (5) with the value T (A) given in Theorem 1.1 shows that when N = 1 the best approximations of T (A) happen when m ∈ {1, 2}. For this reason we have applied the HBM for N ≤ 6 and m ≤ 2 to elucidate which of the approaches is better. In the Table 1 we will compare the percentage of the relative errors e N (m) = 100 The best approximation that we have found corresponds to T 6 (A; 0). Our computers have had problems to get the Gröbner basis needed to fill the gaps of the table. Table 1. Percentage of relative errors e N (m).
The paper is organized as follows. Theorem 1.1 is proved in Section 2. In Section 3 we describe the N-th order HBM adapted to our purposes. Finally, in Section 4 we use this method to demonstrate Theorem 1.2.
(i) We start proving that the solution of (1) with initial conditions x(0) = A, where erf −1 is the inverse of the error function Notice that lim t→± To obtain (6), observe that from system (2) we arrive at the simple differential equation which has separable variables and can be solved by integration. The particular solution that passes by the point (x, y) = (A, 0) is Combining (2) and (7) we obtain again a separable equation. It has the solution (7) we obtain (6), as we wanted to prove.
By using x(t) and y(t) given by (6) and (8), respectively, or using (7), we can draw the phase portrait of (2) which, as we can see in Figure 1.(b), is symmetric with respect to both axes. Notice that its orbits do not cross the y-axis, which is a singular locus for the associated vector field. Moreover, the solutions of (1) are not periodic (see Figure 1.(a)), and the transit time of From (6) we introduce the C 0 -function, defined on the whole R, as 2π}. Hence (i) of the theorem follows. Notice that directly from (1), it is easy to see that this equation can not have C 2 -solutions such that x(t * ) = 0 for some t * ∈ R, because this would imply that lim t→t * ẍ(t) = ∞. (ii) System (3) is Hamiltonian with Hamiltonian function Since ln(x 2 + k 2 ) has a global minimum at 0 and ln(x 2 + k 2 ) tends to infinity when |x| does, system (3) has a global center at the origin. In Figure (3) we can see its phase portrait for some values of k. This figure also illustrates how the periodic orbits of (3) approach to the solutions of system (2). Its period function is Figure 3. Phase portraits of (3) for different values of k.
where h = ln(A 2 + k 2 )/2 is the energy level of the orbit passing through the point (A, 0). Therefore, where we have used the change of variable s = x/A and the symmetry with respect to x. Then, If we prove that and the theorem will follow. Therefore, for completing the proof, it only remains to show that (9) holds. For proving that, take any sequence 1/z n , with z n tending monotonically to infinity, and consider the functions f n (s) = ln . We have that the sequence {f n (s)} n∈N is formed by measurable and positive functions defined on the interval (0, 1). It is not difficult to prove that it is a decreasing sequence. In particular, f n (s) < f 1 (s) for all n > 1. Therefore, if we show that f 1 (s) is integrable, then we can apply the Lebesgue's dominated convergence theorem ( [9]) and (9) will follow. To prove that 1 0 f 1 (s) < ∞ note that for s close to 1, Since this last expression is integrable the result follows by the comparison test for improper integrals.

The Harmonic Balance Method
This section gives a brief description of the HBM applied the second order differential equations F := F (x(t),ẍ(t)) = 0, with F (−u, −v) = F (u, v), and adapted to our interests. Notice that if x(t) is a solution of (10) then x(−t) also is a solution. Suppose that equation (10) has a T -periodic solution x(t) with initial conditions x(0) = A,ẋ(0) = 0 and period T = T (A). If x(t) satisfies x(t) = x(−t) it is clear that its Fourier series has no sinus terms and writes as ∞ k=1 a k cos(kωt), with ∞ k=1 a k = A and ω = 2π T .
As we have seen in previous section, the weak periodic solutions of equation x(t)ẍ(t) + 1 = 0 that we want to approach satisfy the above property. Moreover x(T /4) = 0 andẋ(T /4) does not exist. In any case, if we are searching smooth approximations to this x(t), they should also satisfyẋ(t)ẍ(t) + x(t) ...
x (t) = 0, and henceẋ(T /4) = 0. For this reason, in this work we will search Fourier series in cosinus, not having the even terms cos(2jωt), j ∈ N∪{0}, which do not satisfy this property. This type of a priori simplifications are similar to the ones introduced in [5] for other problems.
3. Find all values a and ω N such that where j N is the value such that (12) consists exactly of N non trivial equations. Notice also that each equation A j (a, ω N , A) = 0 is equivalent to 4. Then the expression (11), with the values of a = a(A) and ω N = ω N (A) obtained in point 3, provide candidates to be approximations of the actual periodic solutions of the initial differential equation. In particular, the functions T N = T N (A) = 2π/ω N give approximations of the periods of the corresponding periodic orbits.
5. Choose, as final approximation, the one associated to the solution that minimizes the norm Remark 3.1. (i) Notice that going from order N to order N + 1 in the method, implies to compute again all the coefficients of the Fourier polynomial, because in general the common Fourier coefficients of x N (t) and x N+1 (t) do not coincide.
(ii) The above set of equations (12) is a system of polynomial equations which usually is not easy to solve. For this reason in many works, see for instance [6,8] and the references therein, only the values of N = 1, 2 are considered. For solving system (12) for N ≥ 3 we use the Gröbner basis approach ( [2]). In general this method is faster that using successive resultants and moreover it does not give spurious solutions.
(iii) As far as we know, the test proposed in point 5 to select the best approach is not commonly used. We propose it following the definition of accuracy of an approximated solution used in [3] and inspired in the classical works [10,11].

Application of the HBM
We start proving a lemma that will allow to reduce our computations to the case A = 1.

Proof. Consider F N = x m+1
Nẍ N + x m N = 0, with x N given in (11). We have to solve the set of N + 1 non-trivial equations with N + 1 unknowns a 1 , a 3 , . . . , a 2N −1 and ω N  Proof of Theorem 1.2. Due to the above lemma, in the application of the N-th order HBM, we can restrict our attention to the case A = 1.
When N = 2, we take an approximation x 2 (t) = a 1 cos(ω 2 t) + a 3 cos(3ω 2 t). The vanishing of the coefficients of 1 and cos(2ω 2 t) in the Fourier series of F 2 provides the nonlinear system By solving it and applying point 5 in the HBM we get that ω 2 = 18/ √ 218. Therefore, as we wanted to prove.
Since all the equations are polynomial, the searching of its solutions can be done by using the Gröbner basis approach, see [2]. Recall that the idea of this method consists in finding a new systems of generators, say G 1 , G 2 , . . . , G ℓ , of the ideal of R[a 1 , a 3 , a 5 , ω 3 ] generated by P, Q, R and S. Hence, solving P = Q = R = S = 0 is equivalent to solve G i = 0, i = 1, . . . , ℓ. In general, choosing the lexicographic ordering in the Gröbner basis approach, we get that the polynomials of the equivalent system have triangular structure with respect to the variables and it can be easily solved. Now, by computing the Gröbner basis of {P, Q, R, S} with respect to the lexicographic ordering [a 1 , a 3 , a 5 , ω 3 ] we obtain a new basis with four polynomials (ℓ = 4), being one of them, Solving G 1 (ω 3 ) = 0 and using again point 5 of our approach to HBM we get that the solution that gives the better approximation is .
The Gröbner basis of {P, Q, R, S, U} with respect to the lexicographic ordering [a 1 , a 3 , a 5 , a 7 , ω 4 ] is a new basis with five polynomials, being one of them an even polynomial in ω 4 of degree 16 with integers coefficients. Solving it we obtain that the best approximation is ω 4 ≈ 1.2453, which gives T 4 (A; 0) ≈ 5.0455 A.
For N = 5 and N = 6 we have done similar computations. In the case N = 5 one of the generators of the Gröbner basis is an even polynomial in ω 5 with integers coefficients and degree 32. When N = 6 the same happens but with a polynomial of degree 64 in ω 6 . Solving the corresponding polynomials we get that ω 5 ≈ 1.2606 and ω 6 ≈ 1.2501, and consequently, T 5 (A; 0) ≈ 4.9843 A, and T 6 (A; 0) ≈ 5.0260 A.
When N = 2, doing similar computations that in item (ii), we arrive at Again, by searching the Gröbner basis of {P, Q, R} with respect to the lexicographic ordering [a 1 , a 3 , ω 2 ] we obtain a new basis with three polynomials, being one of them G 1 (ω 2 ) = 7635411ω 8 2 − 14625556ω 6 2 + 5833600ω 4 2 − 661376ω 2 2 + 13824. Notice that the equation G(ω 2 ) = 0 can be algebraically solved. Nevertheless, for the sake of shortness, we do not give the exact roots. Following again step 5 of our approach we get that the best solution is ω 2 ≈ 1.1915, or equivalently that T 2 (A; 1) ≈ 5.2733A.
Computing the Gröbner basis of {P, Q, R, S} with respect to the lexicographic ordering [a 1 , a 3 , a 5 , ω 3 ] we get that one of the polynomials of the new basis is an even polynomial in ω 3 of degree 26 with integer coefficients. By solving it we obtain that the best approximation is ω 3 ≈ 1.2206, which produces the value T 3 (A; 1) of the statement.
When N = 4 we arrive at five polynomial equations, that we omit. Once more, using the Gröbner basis approach we obtain a polynomial condition in ω 4 of degree 80. Finally, ω 4 ≈ 1.2275 and T 4 (1; A) ≈ 5.1186.
(iv) When m = 2 we have to approach the solutions of F (x(t),ẍ(t)) = x 3 (t)ẍ(t) + x 2 (t) = 0. We do not give the details of the proof because we get our results by using exactly the same type of computations.
Remark 4.2. For each N and m our computations also provide a trigonometric polynomial that approaches the continuous weak periodic solution φ(t) given in the proof of Theorem 1.1.