-2

I wanted to solve this system of differential equations using Python, and I was wondering how to do it. I've tried nothing yet since I'm new to these systems, though I have already solved some individual diff. eqs.. Right now, I'm familiar with solve_ivp and that's practically it.

JoanSGF
  • 101
  • 1
  • 1
    I'd recommend starting by reading a scipy tutorial –  Aug 23 '22 at 14:23
  • Your question reminds me of [The Simpsons](https://www.youtube.com/watch?v=lOTyUfOHgas). Please [do some research](//meta.stackoverflow.com/a/261593/843953) first, this will allow you to formulate a _specific_ question that is [on-topic](/help/on-topic) here. Please take the [tour], and read [ask] and the [question checklist](//meta.stackoverflow.com/q/260648/843953). Once you have done some research and made an attempt, remember to include your [mre] in the question. Welcome to Stack Overflow! – Pranav Hosangadi Aug 23 '22 at 16:15
  • Thanks you both: I finally solved it. I had sworn that I searched on Youtube and here how to solve one of these systems, but I couldn't find any using solve_ivp. I finally found one. I'll post the solution. – JoanSGF Aug 23 '22 at 22:31

1 Answers1

0

Here's the code:

import numpy as np
from numpy import sin,cos,sign,sqrt,pi
from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
from matplotlib.pyplot import figure

g = 9.807
h = 1.8
l = 0.56 * h
R_sc = 14
mu = 0.1
k = 0.3
m = 80
v_0 = 12
phi_0 = np.pi/4

def velphi(t,y):
    dy=np.zeros([3])
    dy[1] = y[2]
    dy[0] = -mu*np.sqrt(g**2 + ((y[0])**2/(R_sc*np.cos(y[1])))**2)-(k/m)*(y[0])**2
    dy[2] = (g/l)*np.sin(y[1])-((y[0])**2/(l*R_sc))*np.sign(y[1])

    return dy

y0 = np.array([v_0, phi_0, 0])
time = np.linspace(0, 10, 10000)
sol = solve_ivp(velphi, (0,10), y0, method='RK45', t_eval=time, dense_output=True, rtol=1e-8, atol=1e-10)
t, pos, vel, phi, omega = sol.t, sol.y[0], sol.y[1], sol.y[2], sol.y[3]

plt.plot(t.T, phi.T,"b",linewidth = 0.65, label = "${\Phi}$")
plt.plot(t.T, omega.T,"g",linewidth = 0.65, label = "${\Omega}$")
s = 'Condicions inicials: ${\Phi}_{o}$ ='+str(np.degrees(phi_0))+'º,  ${\Omega}_{o}$ = 0 rad/s'
plt.title('Pèndol centrífug invertit | '+ s)
plt.xlabel('Temps (s)')
plt.ylabel(u'${\Phi}$ (rad) | ${\Omega}$ (rad / s)')
plt.grid(False)
plt.legend(loc = "upper right")
plt.show()
plt.plot(t.T, vel.T,"r",linewidth = 0.65, label = "v")
s = 'Condicions inicials: ${\Phi}_{o}$ ='+str(np.degrees(phi_0))+'º,  ${\Omega}_{o}$ = 0 rad/s'
plt.title('Pèndol centrífug invertit | '+ s)
plt.xlabel('Temps (s)')
plt.ylabel('v (m/s)')
plt.grid(False)
plt.legend(loc = "upper right")
plt.show()

And here the solutions:

Velocity

Phi and the angular velocity

JoanSGF
  • 101
  • 1