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 omega = r 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) omegapoints.append(r) 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.