Search code examples
pythonfunctionargumentscallableodeint

Passing multiple callable to odeint


I want to write a python script for the following mechanical system :

enter image description here

Where

enter image description here

Which can be rewritten as below :

enter image description here

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.

enter image description here

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, enter image description here

It works, see the result below.

enter image description here

However, when I pass a vector of callable on the RHS it bugs.

enter image description here

enter image description here

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 ?


Solution

  • 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