solving a Non- linear first order differential equation and getting a break in the plot

Question

I am trying to solve the elliptical differential equation using fourth-order runge-kutta method in python.

After execution, I get a very small part of the actual plot that should be obtained and alongside with it an error saying that:

"RuntimeWarning: invalid value encountered in double_scalars"

import numpy as np
import matplotlib.pyplot as plt

#Define constants
g=9.8
L=1.04

#Define the differential Function
def fun(y,x):
    return-(2*(g/L)*(np.cos(y)-np.cos(np.pi/6)))**(1/2)

#Define variable arrays
x=np.zeros(1000)
y=np.zeros(1000)
y[0]=np.pi/6
dx=0.5

#Runge-Kutta Method
for i in range(len(y)-1):
    k1=fun(x[i],y[i])
    k2=fun(x[i]+dx/2, y[i]+dx*k1/2)
    k3=fun(x[i]+dx/2, y[i]+dx*k2/2)
    k4=fun(x[i]+dx, y[i]+dx*k3)
    
    
    y[i+1]=y[i]+dx/6*(k1+2*k2+2*k3+k4)
    x[i+1]=x[i]+dx
    
#print(y)
#print(x)
plt.plot(x,y)
plt.xlabel('Time')
plt.ylabel('Theta')
plt.grid()

And the graph I obtain is something like,

My question is why am I getting the error message? Thanks for helping!

Answer 1

Several points that lead to this behavior. First, you switched the order of the arguments in the ODE function, probably to make it compatible with odeint. Use the tfirst=True optional argument to avoid that and have the independent variable always first.

The actual source of the error is the term

(np.cos(y)-np.cos(np.pi/6)))**(1/2)

remember that in your version y has the value x[i], so that at some point the expression under the root becomes negative.

If you correct the first point, you will probably still encounter the second error as the exact solution moves parabolically towards the fixed point, so that the stages of RK4 are likely to overshoot. One can fix that by providing a sufficiently secured square root function,

def softroot(x): return x/max(1e-12,abs(x))**0.5

#Define the differential Function
def fun(x,y):
    return -(2*(g/L)*softroot(np.cos(y)-np.cos(np.pi/6)))

#Define variable arrays
dx=0.01
x=np.arange(0,1,dx)
y=np.zeros(x.shape)
y[0]=np.pi/6
...

results in a plot

as the solution already starts in the fixed point. Shifting the initial point a little down to y[0]=np.pi/6-1e-8 produces a jump to the fixed point below.