TLDR: I'm solving an ODE that represents the torques applied on a rotating wing that is connected to a motor using scipy.solve_ivp. It seems that the solution doesn't behave as needed, and I can't seem to understand why
Im working on a solution for the equation phi'' = a*torque_{motor}-b*(phi')^2, where a and b are constants and torque_{motor} changes every iteration. More specifically, I wish to solve this ODE in a loop, every time for a differnet time window [t,t+dt],[t+dt,t_2*dt],... where in every such window torque_motor is assigned with a new value (float) - this is like saying that I use different force from the motor every time-step.
A subset of my code is provided as minimal running example:
class Solver():
def __init__(self):
self.torque = 0
def set_torque(self, new):
self.torque = new
def phi_derivatives(self, t, y):
"""
A function that defines the ODE that is to be solved: I * phi_ddot = tau_z - tau_drag.
We think of y as a vector y = [phi,phi_dot]. the ode solves dy/dt = f(y,t)
"""
# for exact values use
a = 1090179.6616082331
b = 271.9406850484702
phi, phi_dot = y[0], y[1]
dy_dt = [phi_dot, a * self.torque - b * (phi_dot ** 2)]
return dy_dt
#
phi_arr = []
phi_dot_arr = []
phi_ddot_arr = []
phi_0 = 0
phi_dot_0 = 0.01
start_t = 0
end_t = 0.05
delta_t = 0.001
s = Solver()
sin = np.sin(2 * np.pi * np.linspace(0, 1, 20))
for torque in sin:
s.set_torque(torque)
sol = solve_ivp(s.phi_derivatives, t_span=(start_t, end_t), y0=[phi_0, phi_dot_0])
phi, phi_dot = sol.y
_, phi_ddot = s.phi_derivatives(0, [phi, phi_dot])
phi_arr.append(phi)
phi_dot_arr.append(phi_dot)
phi_ddot_arr.append(phi_ddot)
phi_0 = phi[-1]
phi_dot_0 = phi_dot[-1]
start_t = end_t
end_t += 0.05
phi_arr = np.concatenate(phi_arr)
phi_dot_arr = np.concatenate(phi_dot_arr)
phi_ddot_arr = np.concatenate(phi_ddot_arr)
notice that for every window I use a torque value that is provided from a sine wave.
Intuitively, I was expecting that the solution phi would oscillate to account for the changing direction of the torque, but what I get is this figure (plotting phi_arr,phi_dot_arr,phi_ddot_arr):
Do you think that this is a matter of bad implementation? also notice the high values of phi' and phi'', which are of order 1e12 and 1e26 respectively