Closed-form answer nonlinear pendulum w/ Jacobi features

The earlier put up seems on the nonlinear pendulum equation and what distinction it makes to the options in the event you linearize the equation.

If the preliminary displacement is sufficiently small, you possibly can merely exchange sin θ with θ. If the preliminary displacement is bigger, you possibly can enhance the accuracy fairly a bit by fixing the linearized equation after which adjusting the interval.

You can even discover a precise answer, however not by way of elementary features; it’s important to use Jacobi elliptic features. These are features considerably analogous to trig features, although it’s not useful to attempt to pin down the analogies. For instance, the Jacobi perform sn is just like the sine perform in some methods however very completely different in others, relying on the vary of arguments.

We begin with the differential equation

θ″(t) + c² sin( θ(t) ) = 0

the place c² = g/L, i.e. the gravitational fixed divided by pendulum size, and preliminary circumstances θ(0) = θ₀ and θ′(0) = 0. We assume −π < θ₀ < π.

Then the answer is

θ(t) = 2 arcsin( a cd(ct | m ) )

the place a = sin(θ₀/2), m = a², and cd is among the 12 Jacobi elliptic features. Notice that cd, like all of the Jacobi features, has an argument and a parameter. Within the equation above the argument is ct and the parameter is m.

The final plot within the earlier put up was deceptive, exhibiting roughly equal components real distinction and error from fixing the differential equation numerically. Right here’s the code that was used to resolve the nonlinear equation.

from scipy.particular import ellipj, ellipk
from numpy import sin, cos, pi, linspace, arcsin
from scipy.combine import solve_ivp

def exact_period(θ):
    return 2*ellipk(sin(θ/2)**2)/pi

def nonlinear_ode(t, z):
    x, y = z
    return [y, -sin(x)]    

theta0 = pi/3
b = 2*pi*exact_period(theta0)
t = linspace(0, 2*b, 2000)

sol = solve_ivp(nonlinear_ode, [0, 2*b], [theta0, 0], t_eval=t)

The answer is contained in sol.y[0].

Let’s examine the numerical answer to the precise answer.

def f(t, c, theta0):
    a = sin(theta0/2)
    m = a**2
    sn, cn, dn, ph = ellipj(c*t, m)
    return 2*arcsin(a*cn/dn)

There are a pair issues to notice in regards to the code. First,SciPy doesn’t implement the cd perform, however it may be computed as cn/dn. Second, the perform ellipj returns 4 features directly as a result of it takes about as a lot time to calculate all 4 because it does to compute one among them.

Here’s a plot of the error in fixing the differential equation.

And right here is the distinction between the precise answer to the nonlinear pendulum equation and the stretched answer to the linear equation.