ODE Pendulum using odeint

Question

I am trying to solve a simple mechanical problem, which is the non linear simple pendulum issue.

For this, I have to use odeint from scipy.

For a reminder, I have to solve the non linar ODE which is: theta'' + w²*sin(theta) = 0

Here is what I have for now, I have to define a function that returns the position of the pendulum, that is to say theta(t):

def ODE(omega,theta0,tf):
    """solution au bout tf avec odeint"""
    
    def F(Y,t,omega):
        """second membre de l'EDO"""
        return np.array([Y[1],-omega**2*np.sin(Y[0])])
    
    Y0 = np.array([theta0,0.])
    Y = odeint(F, Y0, tf, args=(omega,))
    return Y

But when I try a verification with this, to compare to the linear solution:

# test
theta0 = 0.1 
omega = 2.58
N = 15 
T = np.linspace(0,2*np.pi/omega,N)

    
# solution with odeint
ThetaODE = ODE(omega,theta0,T)


# linear solution
ThetaL = theta0*np.cos(omega*T)

# plot
plt.figure(figsize=(15,6))
plt.subplot(1,2,1)
plt.plot(T,ThetaODE,label="ODEint")
plt.plot(T,ThetaL,label="Lineaire")
plt.xlabel('t')
plt.legend()

plt.subplot(1,2,2)
plt.plot(T,ThetaL-ThetaODE)
plt.title('Ecart entre ODEint et la sol Lineaire')
plt.xlabel('t')

I have this error:

ValueError                                Traceback (most recent call last)
<ipython-input-128-d8b2832c3f6c> in <module>
     22 
     23 plt.subplot(1,2,2)
---> 24 plt.plot(T,ThetaL-ThetaODE)
     25 plt.title('Ecart entre ODEint et la sol Lineaire')
     26 plt.xlabel('t')

ValueError: operands could not be broadcast together with shapes (15,) (15,2)

There must be something I understand the wrong way concerning odeint, any idea what I do wrong ? Thank you !

Laurent90 · Answer 1 · 2020-12-02T19:28:14.253

You must understand how the odeint library works. Btw, your should use the library solve_ivp instead: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp

Your system has two variables: theta and its time derivative. Therefore the ODE integration method odeintor solve_ivp will return an array (Yin your function ODE) with two dimensions. One dimension is the number of variables (2), the other is the number of time points that have been computed inbetween (to control that, fere to the documentation).

So here, your plot line should be plt.plot(T,ThetaL-ThetaODE[:,0]).

The advice here is to read the docs and search the web for similar problems. I think you would have found an answer easily. There are tons of tutorials for the integration of ODEs, so please take the time to go through some of them. For example this one has some interesting examples: https://scipy-cookbook.readthedocs.io/items/idx_ordinary_differential_equations.html

This one here shows how to deal with pendulum: https://scientific-python.readthedocs.io/en/latest/notebooks_rst/3_Ordinary_Differential_Equations/04_Exercices/00_Tutorials/02_ODE_tutorial_pendulum.html

Here is the plotting part of your code with the corrections:

# plot the two solution components
plt.figure(figsize=(15,6))
plt.subplot(1,2,1)
plt.plot(T,ThetaODE[:,0],label=r"$\theta$")
plt.plot(T,ThetaODE[:,1],label=r"$\dot{\theta}$")
plt.xlabel('t')
plt.grid()
plt.legend()

# compare the time evolution of the angle with the inear solution
plt.figure(figsize=(15,6))
plt.subplot(1,2,1)
plt.plot(T,ThetaODE[:,0] ,label="ODEint", linestyle='-', marker='+')
plt.plot(T,ThetaL, label="Lineaire", linestyle='--', marker=None)
plt.xlabel('t')
plt.legend()

plt.subplot(1,2,2)
plt.plot(T,ThetaL-ThetaODE[:,0])
plt.title('Ecart entre ODEint et la sol Lineaire')
plt.xlabel('t')

ODE Pendulum using odeint

1 Answers1