将多个可调用对象传递给odeint

Question

我想为以下机械系统编写一个Python脚本：

哪里

可以重写如下：

例如，我们定义以下Python函数来解决此类ODE问题。

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])

上面的脚本有效，cf。如下图.

但是参数是在函数内部给出的，并且给出了强制向量

F = np.array([0, 0, 0, np.cos(3*t), 0, np.cos(2*t)])

。这使得它不可重复使用。

为了使其可重用，我建议使用以下功能。

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

当我使用 RHS 上的数字向量测试脚本时，

成功了，见下面的结果。

但是，当我在 RHS 上传递可调用向量时，它会出现错误。

即使我尝试重塑

u

，也没有改善结果。我尝试将

*callable_args

添加到

ode_dynamic

作为可调用元组的参数，它也不起作用。

我想知道有没有办法将多个可调用对象传递给

odeint

？

Answer 1

经过这里和那里的调整后，我找到了一种将多个可调用对象传递给

odeint

的方法。我对上面的脚本做了一些更改以获得：

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

在

main

即

if __name__==_'__main__':

中，右侧给出了一组

lambda

F1 = lambda t : np.cos(3*t)
F2 = lambda t : 0
F3 = lambda t : 2*np.cos(3*t)

并传递给

odeint

，如下所示：

X_sol = odeint(ode_dynamic, y0, t, args=(M, C, K, F1, F2, F3), tfirst=True)

这就是如何成功地将多个时间相关函数传递给

odeint

，即使右侧是一组数字，它也必须定义为

callable

F1 = lambda t : 6.0
F2 = lambda t : 0.0
F3 = lambda t : -2.0
F4 = lambda t : np.pi