Attempt to solve the nonlinear pendulum 2nd order differential equation using 4th order Runge-Kutta method, not getting expected result

Question

I am new to Python, so I have restricted the process to matplotlib only, not going into NumPy.

I followed the book "Scientific Computing in Python" by Abhijit Kar Gupta, to write the Python script to solve a 2nd-order differential equation using the 4th order Runge-Kutta method.

In the book, the example the author gave was a simple one. But for my problem, the nonlinear pendulum (where small angle approximation sin(x)~x is inaccurate), I tried to modify the functions as I expected would be necessary.

I am trying to plot the variation of the angle (x) vs. time (t). But I don't think I am getting the correct output. The nature of a pendulum released at an angle of x = 170° is given here: https://en.wikipedia.org/wiki/Pendulum_(mechanics)#/media/File:Pendulum_170deg.gif.

Though I'm finding it for x = 160° = 8.0*pi/9.0 radians, still the plot should look somewhat like a sinusoidal oscillation, with a non-uniform period. But all I'm getting is something like a suddenly stopped pendulum.

Please tell me where I am doing it wrong. Also, what should I do if, rather than looking at the plot, I want to get a list-like output, like a table with two columns for time and angle? Also also, what should be done so that we shall get the angle for a particular time (that is, for any external input)?

Here is the code I wrote:

Python 3.12.0 (tags/v3.12.0:0fb18b0, Oct  2 2023, 13:03:39) [MSC v.1935 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license()" for more information.

# Solving a non-linear pendulum equation of motion using RK4 method:
# Equation of Motion: d^x/dt^2 + k.sin(x)=0, where 'x' is angle, k=g/l. We take g=10m/s^2, l=10m, so k=1(s^{-2}).
# We modify the single 2nd order differential equation as a system of linear equations, consisting of two 1st order differential equation: (1) dx/dt=y=f1(x,y,t), (2) dy/dt=-k.sin(x)=f2(x,y,t).

from math import sin, pi

# Defining two functions-
def f1(t,x,y): return y

def f2(t,x,y): return -k*sin(x)

k=1.0                      #parameter
t,x,y=0,8.0*pi/9.0,0             #initial values (t: second, x: radian, y: (initial angular velocity) radian/second)
h=0.01                   #increment in t
T,X,Y=[t],[x],[y]          #lists to store data

# Loop starts with RK4 steps:
for i in range(200):
    a1=h*f1(t,x,y)                      
    b1=h*f2(t,x,y)                      
    a2=h*f1(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    b2=h*f2(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    a3=h*f1(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    b3=h*f2(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    a4=h*f1(t+h,x+a3,y+b3)              
    b4=h*f2(t+h,x+a3,y+b3)              
    x=x+(1/6)*(a1+2*a2+2*a3+a4)         # apprx. value of x1 for t+h
    y=y+(1/6)*(b1+2*b2+2*b3+b4)         # apprx. value of y1 for t+h
    t=t+h                               # current value of independent variable t
    T.append(t)
    X.append(x)
    Y.append(y)

    
# for plotting
import matplotlib.pyplot as plt
plt.plot(T,X,lw=3)
[<matplotlib.lines.Line2D object at 0x00000190D0CD5910>]
plt.xlim(0,50)    # xlim is just time range
(0.0, 50.0)
plt.ylim(-pi,pi)  # ylim is the range of angle (in radians)
(-3.141592653589793, 3.141592653589793)
plt.grid()
plt.show()

And the plot I got:

If I zoom the blue line:

This is in no way close to a general pendulum. Why is the solution stopping abruptly?

Answer 1

As have been pointed out in the comments, the problem is that you chose only 20 samples for plotting. You can parameterize the number of samples as a function of the time, and thus get a graph where t ranges from 0 to 50, as you intended:

# Solving a non-linear pendulum equation of motion using RK4 method:
# Equation of Motion: d^x/dt^2 + k.sin(x)=0, where 'x' is angle, k=g/l. We take g=10m/s^2, l=10m, so k=1(s^{-2}).
# We modify the single 2nd order differential equation as a system of linear equations, consisting of two 1st order differential equation: (1) dx/dt=y=f1(x,y,t), (2) dy/dt=-k.sin(x)=f2(x,y,t).

from math import sin, pi

# Defining two functions-
def f1(t,x,y): return y

def f2(t,x,y): return -k*sin(x)

k=1.0                      #parameter
t,x,y=0,8.0*pi/9.0,0             #initial values (t: second, x: radian, y: (initial angular velocity) radian/second)
h=0.01                   #increment in t
T,X,Y=[t],[x],[y]          #lists to store data
tmax=50                  #max value of t (duration of the plot in units of time)
n=int(tmax/h)            #number of samples

# Loop starts with RK4 steps:
for i in range(n):
    a1=h*f1(t,x,y)                      
    b1=h*f2(t,x,y)                      
    a2=h*f1(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    b2=h*f2(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    a3=h*f1(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    b3=h*f2(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    a4=h*f1(t+h,x+a3,y+b3)              
    b4=h*f2(t+h,x+a3,y+b3)              
    x=x+(1/6)*(a1+2*a2+2*a3+a4)         # apprx. value of x1 for t+h
    y=y+(1/6)*(b1+2*b2+2*b3+b4)         # apprx. value of y1 for t+h
    t=t+h                               # current value of independent variable t
    T.append(t)
    X.append(x)
    Y.append(y)

    
# for plotting
import matplotlib.pyplot as plt
plt.plot(T,X,lw=3)
plt.xlim(0,tmax)    # xlim is just time range
plt.ylim(-pi,pi)  # ylim is the range of angle (in radians)
plt.grid()
plt.show()

Also, what should I do if, rather than looking at the plot, I want to get a list-like output, like a table with two columns for time and angle?

You already have the values for time and angle in variables T and X, respectively. If you don't want to use numpy, you could combine those arrays in a matrix-like structure like so:

TX = [[t,x] for t,x in zip(T,X)]

As to your following question:

Also also, what should be done so that we shall get the angle for a particular time (that is, for any external input)?

For an input that happens to be in your array T, you only have to get the corresponding element in array X. For an arbitrary value of time, you could do a linear interpolation between the closest samples. This could be more conveniently done using a function:

def get_angle(t, T, X, h):
    # t: Any particular time (must be in the range used for T)
    # T: The time array you calculated earlier
    # X: The angle array you calculated earlier
    # h: Your time step (increment in t)

    i = t/h   # The corresponding sample index (not necessarily an integer)
    i0 = floor(i)
    i1 = ceil(i)
    x = (X[i1]-X[i0])*(i-i0)/(i1-i0) + X[i0]    # Linear interpolation between the two closest samples
    return x