我试图用sympy解决一个复杂的微分方程组 . 我使用sympy来快速计算时间导数，之后我有两个包含导数本身的导数方程 . 方程式不是线性的，不适 Contract 意识别的方程式，因此它会抛出一个未实现的错误 . 有没有更简单的方法来解决这些方程（甚至数字），并获得各自的运动定律（跨时间的值）？这可能非常低效，所以如果有人知道一个更有效的过程，我很满意，我主要是开始使用sympy，因为计算衍生物相当快（否则我会花很多时间浪费纸张） .
import sympy as sym
a,b = sym.S(['a','b'])
S1,S2, alpha, r, c1 ,c2 = sym.symbols('S1, S2, alpha, r, c1, c2', negative=False)
t = sym.var('t')
x1 = sym.Function("x1")(t)
x2 = sym.Function("x2")(t)
lam = sym.Function('lam')(t)
gam = sym.Function('gam')(t)
p = (1/3)*a*(alpha*x1 + (1-alpha)*x2)**3 + (1/2)*b*(alpha*x1 + (1-alpha)*x2)
lagrangian = p - c1/2*alpha*x1 - c2/2*(1-alpha)*x2 + lam*(S1 - 0.5*x1**2) + gam*(S2 - 0.5*x2**2)
FOC_x1 = sym.diff(lagrangian,x1)
FOC_x2 = sym.diff(lagrangian,x2)
x1_lam = sym.solve(FOC_x1,x1)[0]
lam_x1 = sym.solve(FOC_x1,lam)[0]
x2_gam = sym.solve(FOC_x2,x2)[0]
gam_x2 = sym.solve(FOC_x2,gam)[0]
dx1_dt = sym.diff(x1_lam,t).subs(lam,lam_x1)
dx2_dt = sym.diff(x2_gam,t).subs(gam,gam_x2)
x1dot = sym.Derivative(x1,t)
x2dot = sym.Derivative(x2,t)
eq = [sym.Eq(x1dot,dx1_dt),sym.Eq(x2dot,dx2_dt)]
sym.dsolve(eq)
为了澄清：我的状态变量分别是S1和S2，我的控制变量分别是x1和x2，最后我的costate变量分别是lam和gam .
1 回答
我不太能确定你的因变量是什么，
x1和x2？对于表示为SymPy表达式的一阶常微分方程（系统）的数值积分，您可能有兴趣使用pyodesys中的SymbolicSys（示例找到here） .