# Solving nonlinear pendulum using Runge-Kutta 4 for smaller steps

I am trying to solve nonlinear pendulum using 4th order Runge-Kutta method for limits between a=0.0 to b=110 seconds and simulated the results to observe the pendulum movement. But when I increase the number of steps N in RK4 method, the pendulum isn't decaying, it is oscillating forever, but for smaller value of N, I can observe the decaying effects.

from math import sin, pi
import numpy as np
import matplotlib.pyplot as plt

# constants
G = 9.51        # acceleration due to gravity
L = 0.1         # length of the pendulum
T = 0.5         # time period of the pendulum
ANGLE = 90.0    # initial angle of the pendulum

def f(r, t):
"""Vectorized simultaneous ODEs functions."""
theta = r[0]
omega = r[1]
ftheta = omega
fomega = - (G / L) * sin(theta)

return np.array([ftheta, fomega], dtype=float)

a = 0.0
b = 110.0
N = 2000
h = (b - a) / N

tpoints = np.arange(a, b, h)
thetapoints = []
omegapoints = []

theta0 = (pi / 180) * ANGLE
omega0 = (2 * pi) / T
r = np.array([theta0, omega0], dtype=float)
for t in tpoints:
thetapoints.append(r[0])
omegapoints.append(r[1])
k1 = h * f(r, t)
k2 = h * f(r + 0.5 * k1, t + 0.5 * h)
k3 = h * f(r + 0.5 * k2, t + 0.5 * h)
k4 = h * f(r + k3, t + h)
r += (k1 + 2 * k2 + 2 * k3 + k4) / 6

ypoints = -L * np.cos(thetapoints)
xpoints = L * np.sin(thetapoints)

plt.plot(tpoints, omegapoints, 'k', lw=0.8)
plt.xlabel('t')
plt.ylabel('omega')
plt.savefig('2000_steps.jpeg', dpi=300, format='jpeg')
plt.show()



Output of omega vs time:

N = 1000

N = 2000

I strongly believe this is due to the approximation

$$\sin \theta \approx \theta$$

for smaller values of $$\theta$$ in the ordinary differential equation.

• in which way do you think is related the "small-angles" approximation with the RK scheme? – VoB Apr 30 '20 at 10:38
• @VoB I think for smaller values of h << 1, function are computed in RK4 which results in smaller value of theta leading to $\sin \theta \approx \theta$. Theta plots are similar to omega plots. – 147875 Apr 30 '20 at 10:43
• Yes, this is one point. The other one I think is related with the fact that RK is not A-stable (and so neither L-stable). So I suggest you to try with an implicit method like backward euler: you should see the correct behaviour – VoB Apr 30 '20 at 11:37