Consider the simple-looking second-order differential equation

d2g dt2 +gn=1

where the exponent n can take any positive real value. This equation appropriately scaled describes the evolution of the square of the radial variable of a body moving under the influence of a power potential.

The equation is in general hopelessly nonlinear and as such not amenable to exact solution. It turns out, however, that a very good approximate solution can be found for continuous variation of the exponent. Initial conditions for the solution will have a vanishing first derivative at the origin, which restricts the solution to functions even in the variable of evolution. This corresponds to a body beginning movement from one turning point.

Numerical solutions with these initial conditions look very similar as a function of the exponent, as can be readily seen in this interactive graphic:

The initial condition c0 is limited by the value n of the exponent in a manner to be described later.

For all values of the exponent there is one trivial exact solution: g(t)=1  . This solution corresponds to a circular orbit at the maximum angular momentum allowed for bound orbits.

For positive values of the exponent there are two simple exact solutions. For n = 1 the solution is a simple sinusoidal,

g(t) =1+(c0 -1)cost , n=1

as can be verified by differentiation. The solution is symmetric about the trivial solution, with extrema at 1+(c0 -1) and 1-(c0 -1)  . Using a half-angle identity this can also be written

g(t) =2-c0 +2(c0 -1) cos2t2 , n=1

For n = 2 the exact solution is given by elliptic functions. Since any elliptic function can be expressed in terms of other elliptic functions with the same periods, choose a representation of the solution that relates easily to the form just given:

g(t) =c1 +(c0 -c1) cn2 (ωt,k) , n=2

The solution is no longer symmetric about the trivial solution, necessitating two constants for the extrema. The second of the two extrema occurs at quarter periods of the Jacobi elliptic function,

ωt1 =K(k) cn[K(k), k] =0

where the elliptic function becomes zero. The second derivative of the squared Jacobi elliptic function is

d2 dt2 cn2 (ωt,k) =2ω2[ sn2 (ωt,k) dn2 (ωt,k) -cn2 (ωt,k) dn2 (ωt,k) +k2sn2 (ωt,k) cn2 (ωt,k) ] d2 dt2 cn2 (ωt,k) =2ω2[ 1-k2 +2(2k2 -1)cn2 (ωt,k) -3k2 cn4 (ωt,k) ]

Substituting this in the differential equation leads to a set of equations among the four constants,

2ω2 (1 -k2) (c0 -c1) +c12 =1 2ω2 (2k2 -1) +c1 =0 6ω2 k2 +(c0 -c1) =0

Using the last two equations to eliminate the extrema from the first equation gives

4ω4 (k4 -k2 +1) =1

while adding the last two equations together gives the relation

2ω2 (k2 +1) =c0

Eliminating the frequency between these two relations, the elliptic modulus is a solution of

(c02 -1)k4 -(c02 +2)k2 +c02-1 =0 k2 =c02 +2±3 c04 -c02 2(c02 -1) =[3c0 ±4 -c02 ]2 4(c02 -1)

The square of the frequency is

ω2 =c0 2(k2 +1) =c02 -1 3[ 3c0 ±4 -c02 ]

and the second extremum is

c1 =c0 -6ω2 k2 =c0 -32 [3c0 ±4 -c02 ]

Since both extrema are equal to unity for the trivial exact solution, the correct sign for the radical is the lower one. The final solution, with a bit of simplification, is thus

c1 =c02 +32 4 -c02 3 ω=12 c0 +4 -c02 3 k=3 c0 -4 -c02 2c02 -1 =2 c02 -1 3c0 +4 -c02

The frequency and elliptic modulus have been rationalized to show that they both remain finite for the trivial solution, with the latter going to zero.

For an solution analytic in the exponent n, the constants here would include a functional dependence on that value. Since it is a property of the Jacobi elliptic function that cn(x,k) cosx for zero modulus, the exact solution for n = 2 has a simple relationship to that for n = 1 . If the functional dependence on exponent of the former solution were to allow k → 0 for arbitrary c0 as n → 1 , then the exact solution for n = 1 would then be recovered with the additional limits

c1 2-c0 and ω12 , n1

While such a solution analytic in the exponent n is not forthcoming, the relationship between these two exact solutions suggests that one of them can be used as an approximate solution for all values of the exponent. Since the the squared Jacobi function already displays an asymmetry about the trivial solution that the circular function does not, the former is a more appropriate choice for the approximation.

Of the three constants needed apart from the initial condition to describe the approximate solution, two can be stated exactly. First multiply the original differential equation by the first derivative of the dependent variable and integrate once:

dg dt [d2g dt2 +gn=1] 12 (dg dt )2 +g n+1 n+1 =g-l

The constant of integration l is proportional to the square of angular momentum in physical applications, hence the letter. At the extrema of the motion the first derivative goes to zero so that

c0 -c0 n+1 n+1 =l =c1 -c1 n+1 n+1

The second extremum is thus a solution to the non-polynomial equation

c1 n+1 n+1 -c1+l =0

With the scaling c1 =(n+1 )1/n c this becomes

cn+1 -c +l(n+1 )1/n =0

which is close to the functional form given at the beginning of a presentation on the real roots of non-polynomials. Accounting for minor differences, the solution found for the root in the vicinity of zero is

czero = i=0 Γ[ni +i+1] Γ[ni +2] 1i! [l(n+1 )1/n ]ni +1

and the second extremum in terms of the first is

c1 =(n+1 )1/n i=0 Γ[ni +i+1] Γ[ni +2] 1i! [ (n+1) c0 -c0 n+1 (n+1 )1 +1/n ]ni +1

It was demonstrated in the presentation that this series converges for its entire domain of application, so that in principle it can be taken as an exact result. In practice it is simpler and more accurate to evaluate the second extremum as a numerical root of the non-polynomial equation.

The real periodic solutions of the differential equation reach a limit when the second extremum is zero, since an arbitrary power of a negative functional value would necessitate a complex solution. Inspection of the series above indicates that the corresponding value for the first extremum occurs when each term is zero, or

(c0 )max =(n+1 )1/n

which is the limitation mentioned at the outset. This limit corresponds to an orbit with zero angular momentum.

The bracketing values of this limit can be found by taking a logarithm and applying l’Hôpital’s rule,

ddn ddn log(n+1) n =1n+1

so that reversing the logarithm with an exponential gives

limn0 (c0 )max =e limn (c0 )max =1

For values of the initial condition larger than this limit the constant l becomes negative, and being proportional to the square of angular momentum implies complex values for the latter. This potentially interesting feature of the equation will be left for a future presentation.

Formally integrating the differential equation provides an expression for the independent variable:

t=12 dg g- gn+1 n+1 -l

Integrating between the two extrema gives the time between, which as stated above is a quarter period of the Jacobi elliptic function. Replacing g with a dummy variable for clarity, the frequency of the Jacobi elliptic function in terms of the initial condition and the elliptic modulus is

ω=K(k) t1 =2 K(k) [ c1c0 dc c- cn+1 n+1 -l ]1

The order of endpoints in the integral keeps the frequency positive. This is an exact expression for a given elliptic modulus k and ensures that the approximation always matches an exact solution of the differential equation at extrema. This can be seen by comparing the approximation so far with the numerical integration given at the outset:

Even though the denominator of the integral in the frequency is zero at both extrema, the integral can be numerically evaluated quite handily using tanh-sinh quadrature, apart from when c0 = c1 = 1 . In this limit integration is problematic and can be replaced with an exact result. If the extrema are written as c ≈ 1 ± δ , then applying the binomial series to the first two terms of the denominator gives

c-c n+1 n+1 1±δ -1n+1 [1±(n +1)δ +(n+1) n2 δ2] c-c n+1 n+1 nn+1 -n2δ2 =nn+1 -n2 (c-1 )2

Leaving the third term in the denominator initially untouched, the integral can be evaluated as an inverse sine. The constant l, being defined by the first two terms of the denominator at the extrema, is then replaced using the same expansion for the final result:

1-δ 1+δ dc c-c n+1 n+1 -l 1-δ 1+δ dcn n+1 -l -n2 (c-1 )2 =δ δ du nn+1 -l -n2 u2 =2n sin1[ un2 nn+1 -l] |δ δ =2n sin1[ uδ] |δ δ =2n π

The approximation above in red is to the eye already a pretty good match for the numerical solution for arbitrary elliptic modulus. The discrepancy can be seen better by visualizing the difference between the numerical and approximate solutions:

Varying the elliptic modulus k shows that the approximation, while matching the numerical solution at extrema, is generally either slightly above or slightly below the latter between extrema. Careful variation with cursor keys will show that there is a small range of elliptic modulus for which the approximation appears to be equally above and below the numerical solution between extrema.

To make this more exact, define a simple deviation by

σ= evaluation points [gnumerical (t) -gapproximate (t) ]2

where the summation is taken over some number of points equally spaced between extrema. The value of this deviation will of course depend the number of points used, but its value is less important than its general behavior. Choosing the number of points equal to ten for calculational expediency, deviation as a function of elliptic modulus looks like this:

Clearly there is a value of elliptic modulus for which the deviation assumes a minimum, at least for n ≥ 1 . The location of the minimum is insensitive to the number of points chosen.

A natural best choice for the elliptic modulus would thus appear to be the value that minimizes the deviation for a given exponent and initial condition. Estimating the behavior of these values as a function of the two other parameters will complete the approximate solution of the differential equation.

The behavior of the values is unclear for n < 1 because there is only an endpoint minimum at k = 0 . Since this single value is unlikely to give an accurate description of the solution over the whole range 0 < n < 1 , the remainder of the presentation will be restricted to n ≥ 1 .

Directly following this paragraph is a hidden interactive graphic ( ) that evaluates locations of minima for given n as a function of c0. Since the evaluations take some time, they have already been processed with this graphic and the resulting data embedded in the page. This allows subsequent visualizations of the data to be immediately interactive. The embedded data is explicitly constrained to be identically zero for both n = 1 and c0 = 1 : the former condition matches the behavior of the exact solution for n = 1 , while the second is consistent with the behavior of the exact solution for n = 2 .

With the minima in deviation evaluated and embedded, visualizing their behavior is straightforward. In three dimensions this looks like

but is perhaps more usefully visualized in two dimensions as a function of exponent for a given initial condition,

or as a function of initial condition for a given exponent,

along with the behavior in exponent for the maximum allowed initial condition c0 =(n+1 )1/n  :

For most of their domains, the first two graphs are both reminiscent of an inverse secant, with the third mostly definitely so. All of these graphs go to zero for either n = 1 or c0 = 1 , and since sec1 (1) =0 the desired behavior is achieved by having an argument in (n - 1) and (c0 - 1) offset by unity. The graphs also appear to have asymptotes to unity for large values, and this is achieved by dividing by the asymptotic value of the inverse secant. The behavior of the elliptic modulus that minimizes the deviation from the numerical solution is then of the form

k(n,c0) 2π sec1 [f1( n-1) f2( c0-1) +1]

where f1(0) =f2(0) =0  . Since the solution to this point is approximate and not exact, what is needed now are choices for these two functions that do a decent job of matching the available data. The functions should be as simple as possible, without worrying about finding a precise match to the data.

Begin with the third graph, the behavior in exponent for the maximum allowed initial condition, since it has fewer variables. One estimate arises from temporarily removing the dependence on c0 by setting f2=1 and letting the other function be a simple linear: f1 =a(n-1)  . Visualizing a variation of the single constant

indicates that there is a value between one and two where the two lines coincide quite nicely. A convenient way to set this value is from the exact solution for n = 2 , for which the limit value c0=3 and so

k=12 =2π sec1 (a+1) a=sec( π22 ) -1 1.252

The remaining behavior for the argument of the inverse secant can be approximated by taking

f2=[ c0-1 (n+1 )1/n -1 ]b

The need for the exponent b can be seen by visualizing the dependence of elliptic modulus k as a function of exponent n along with the numerical data,

or as a function of initial condition for a given exponent along with the numerical data:

There again appears to be a value between one and two for which the estimate matches the numerical data fairly well, albeit not nearly as exactly as that for the behavior in exponent for the maximum allowed initial condition. Since there is no compelling way to select this value by comparison to known exact solutions, a value of b = 1.75 will be chosen. To the same level of approximation, let a = 1.25 .

With these approximate values for the parameters, the three-dimensional behavior in comparison with the numerical data is

which is a pretty good approximation to the surface, considering the complexity of the actual problem.

In summary, the differential equation

d2g dt2 +gn=1

has a decent approximate solution of the form

g(t) c1 +(c0 -c1) cn2 (ωt,k)

with the three dependent parameters given by

c1 : c1 n+1 n+1 -c1+l =0 k2π sec1 [1.25( n-1) [ c0-1 (n+1 )1/n -1 ]1.75 +1] ω=2 K(k) [ c1c0 dc c- cn+1 n+1 -l ]1 , l=c0 -c0 n+1 n+1

In order for the solution to remain real, the initial condition lies between unity and

(c0 )max =(n+1 )1/n

Keep in mind that the behavior of the elliptic modulus was constructed to minimize the deviation of the approximation from the numerical solution. When applying this approximate solution to the orbits determined by power potentials with arbitrary real exponents, it will turn out that another choice is more appropriate.

Uploaded 2016.12.27 — Updated 2017.09.30