Python求解含二阶导数项的2自由度气动弹性ODE系统问题

你的问题核心在于Matlab和Python对变量/数组元素赋值的处理逻辑不同：Matlab中数组元素的赋值可以互相引用（延迟计算），但Python要求变量必须先定义才能被引用，所以直接照搬Matlab的写法会出现UnboundLocalError。

解决方法是联立求解二阶导数的线性方程组，把两个相互依赖的hdd和alphadd当作二元一次方程组的未知数，直接计算解析解。

步骤1：整理原ODE为线性方程组 #

原方程：

m*h'' + S*alpha'' = L(t) - k_h*h
S*h'' + I*alpha'' = M(t) - k_alpha*alpha

令hdd = h''，alphadd = alpha''，写成线性方程组形式：

[ m   S ] [hdd   ]   = [ L - k_h*h     ]
[ S   I ] [alphadd]     [ M - k_alpha*alpha ]

步骤2：手动计算解析解（无需额外库） #

利用克莱姆法则求解：

• 系数矩阵行列式：det_A = m*I - S**2
• 计算hdd和alphadd：

  hdd = (I*(L - k_h*h) - S*(M - k_alpha*alpha)) / det_A
  alphadd = (m*(M - k_alpha*alpha) - S*(L - k_h*h)) / det_A
  

步骤3：Python代码实现 #

手动计算版本 #

def ode(t, y, m, I, S, k_h, k_alpha, L, M):
    h, hd, alpha, alphad = y
    
    # 计算行列式
    det_A = m * I - S ** 2
    # 计算右边项
    b1 = L - k_h * h
    b2 = M - k_alpha * alpha
    
    # 求解二阶导数
    hdd = (I * b1 - S * b2) / det_A
    alphadd = (m * b2 - S * b1) / det_A
    
    return [hd, hdd, alphad, alphadd]

使用numpy求解（更适合复杂系统） #

如果后续扩展自由度，用numpy的线性方程组求解更简洁：

import numpy as np

def ode(t, y, m, I, S, k_h, k_alpha, L, M):
    h, hd, alpha, alphad = y
    
    # 定义系数矩阵和右边向量
    A = np.array([[m, S], [S, I]])
    b = np.array([L - k_h*h, M - k_alpha*alpha])
    
    # 求解方程组
    hdd, alphadd = np.linalg.solve(A, b)
    
    return [hd, hdd, alphad, alphadd]

为什么Matlab可以运行？ #

Matlab中dYdt是一个预先分配的数组（即使未显式初始化，也会在第一次赋值时创建），赋值dYdt(2) = ... dYdt(4)...时，Matlab会保留这个表达式的引用，等dYdt(4)赋值完成后，整个数组的元素会被正确计算。而Python是逐行执行的变量赋值，hdd = ... alphadd ...时alphadd还未定义，因此报错。

内容的提问来源于stack exchange，提问作者Gabriel Araujo