0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ in which way do you think is related the "small-angles" approximation with the RK scheme? $\endgroup$ – VoB Apr 30 at 10:38
  • $\begingroup$ @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. $\endgroup$ – 147875 Apr 30 at 10:43
  • $\begingroup$ 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 $\endgroup$ – VoB Apr 30 at 11:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.