I want to write a python script for the following mechanical system :
Where
Which can be rewritten as below :
For example, we define the following Python function to solve this kind of ODE problem.
def ode_dynamic(t, Y):
# system parameters
masse1, masse2, masse3 = 550, 450, 850
k_raideur1, k_raideur2, k_raideur3, k_raideur4 = 450, 70, 350, 510
c_amort1, c_amort2, c_amort3, c_amort4 = 10, 20, 40, 5
M = np.diag([masse1, masse2, masse3])
K = np.array([[k_raideur1 + k_raideur2, -k_raideur2, 0],
[-k_raideur2, k_raideur2 + k_raideur3, -k_raideur3],
[0, -k_raideur3, k_raideur3 + k_raideur4]])
C = np.array([[c_amort1 + c_amort2, -c_amort2, 0],
[-c_amort2, c_amort2 + c_amort3, -c_amort3],
[0, -c_amort3, c_amort3 + c_amort4]])
y1, y2, y3, y4, y5, y6 = Y
A = -np.matmul(np.linalg.inv(M), K)
B = -np.matmul(np.linalg.inv(M), C)
D = np.block([[np.zeros([len(M), len(M)]), np.eye(len(M))],[A, B]])
dy1 = y4
dy2 = y5
dy3 = y6
dy4 = D[3, 0]*y1 + D[3, 1]*y2 + D[3,2]*y3 + D[3,3]*y4 + D[3,4]*y5 + D[5,4]*y6 + np.cos(3*t)
dy5 = D[4, 0]*y1 + D[4, 1]*y2 + D[4,2]*y3 + D[4,3]*y4 + D[4,4]*y5 + D[5,4]*y6
dy6 = D[5, 0]*y1 + D[5, 1]*y2 + D[5,2]*y3 + D[5,3]*y4 + D[5,4]*y5 + D[5,4]*y6 + np.cos(2*t)
return np.array([dy1, dy2, dy3, dy4, dy5, dy6])
The above script works, cf. figure below.
But the parameters were given inside the function and the forcing vector is given F = np.array([0, 0, 0, np.cos(3*t), 0, np.cos(2*t)]). Which makes it not reusable.
To make it reusable, I propose the function below.
def ode_dynamic(t, y, M, C, K, F):
A = -np.matmul(np.linalg.inv(M), K)
B = -np.matmul(np.linalg.inv(M), C)
D = np.block([[np.zeros([len(M), len(M)]), np.eye(len(M))],[A, B]])
'''
Setting the right hand side
'''
F_zeros = np.zeros(len(M))
U = np.block([F_zeros, np.matmul(np.linalg.inv(M), F)])
# returning
return np.matmul(D, y) + U
When I test the script with a vector of numbers on the RHS,
It works, see the result below.
However, when I pass a vector of callable on the RHS it bugs.
Even when I tried to reshape the u, it didn’t improve the result. I tried to add *callable_args to ode_dynamic as parameters of a tuple of callable, it didn’t work either.
I am wondering is there a way to pass multiple callable to odeint ?
After tweaking here and there, I found a way to pass multiple callable to odeint. I made minor changes in the script above to get the script below :
def ode_dynamic(t, y, M, C, K, *callable_args):
A = -np.matmul(np.linalg.inv(M), K)
B = -np.matmul(np.linalg.inv(M), C)
D = np.block([[np.zeros([len(M), len(M)]), np.eye(len(M))],[A, B]])
# Unpack callables
F_eval = np.array([f(t) for f in callable_args])
'''
Setting the second member of the function
'''
u_zeros = np.zeros(len(M))
U = np.block([u_zeros, np.matmul(np.linalg.inv(M), F_eval)])
return np.matmul(D, y) + U
And in the main i.e. if __name__==_'__main__': the right-hand side is given as a set of lambda
F1 = lambda t : np.cos(3*t)
F2 = lambda t : 0
F3 = lambda t : 2*np.cos(3*t)
And passed to the odeint like the following :
X_sol = odeint(ode_dynamic, y0, t, args=(M, C, K, F1, F2, F3), tfirst=True)
This is how I succeeded in passing multiple time dependent function to odeint even if the right-hand side is a set of numbers it must be defined as callable
F1 = lambda t : 6.0
F2 = lambda t : 0.0
F3 = lambda t : -2.0
F4 = lambda t : np.pi