从渐近性出发求解二阶微分方程组 [英] Solving a system of 2nd order differential equations from sympy

问题描述

我正在做一个多自由度动力学问题，使用的是二阶拉格朗日方程。我用渐近解出了运动方程。现在这些方程在计算导数后变得相当长，尽管符号简化不能进一步简化它。我的问题实际上是如何从这里解决这个三个二阶颂歌的系统。我不知道如何转换这些公式，这样它们才能与scipy.odeint()一起使用。我想到了替换，但有很多符号。所以我在搜索phi0，phi1和phi2，以及它们的一阶和二阶导数。初始条件都是phi[0]=0，所有dphi[0]=0。 我希望有一种方法可以解决这个问题，而不是从头开始。先谢谢你。

def derivativeLagranga(Lagrange,n):
"""left side of lagrange"""
f0 = sym.Function('f0')(t)
f1 = sym.Function('f1')(t)
f2 = sym.Function('f2')(t)
f3 = sym.Function('f3')(t)
L_i = []
L_it = []
L_j =[]
L_leva = []
x=0
y=0
for i in range(0,n-1):
    x = Lagrange.diff(kot[i].diff(t))
    L_i.append(x)
for i in range(0,n-1):
    x = L_i[i].diff(t)
    x = x.replace(sym.sin(kot[i]),kot[i])
    L_it.append(x)
for i in range(0,n-1):
    x = L.diff(kot[i])
    L_j.append(x)

for i in range(0,n-1):
    x = L_it[i]+L_j[i]
    L_leva.append(x)



return L_left

left_side_L = derivativeLagranga(Lagrange, n)

f0 = sym.simplify(left_side_L[0].subs(values))
f1 = leva_stran_L[1].subs(values)
f2 = leva_stran_L[2].subs(values)

f0

所以我的f0是其中的一个方程，我不能复制输出，所以我会发布一张图片。

54.51345(−(𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+2.0𝜑0(𝑡)𝑑𝑑𝑡𝜑0(𝑡))(0.5(𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))cos(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))−1.0cos(𝜑0(𝑡))𝑑𝑑𝑡𝜑0(𝑡)−1.0cos(𝜑1(𝑡))𝑑𝑑𝑡𝜑1(𝑡)+1.0cos(𝜑2(𝑡))𝑑𝑑𝑡𝜑2(𝑡))+54.51345((𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))cos(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+cos(𝜑0(𝑡))𝑑𝑑𝑡𝜑0(𝑡))(0.5(𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+0.5𝜑0(𝑡)𝑑𝑑𝑡𝜑0(𝑡)+0.5sin(𝜑1(𝑡))𝑑𝑑𝑡𝜑1(𝑡)+1.0sin(𝜑2(𝑡))𝑑𝑑𝑡𝜑2(𝑡))+54.51345(2.0(𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))cos(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+cos(𝜑0(𝑡))𝑑𝑑𝑡𝜑0(𝑡))(1.0(𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+0.5𝜑0(𝑡)𝑑𝑑𝑡𝜑0(𝑡)+0.5sin(𝜑1(𝑡))𝑑𝑑𝑡𝜑1(𝑡)+1.0sin(𝜑2(𝑡))𝑑𝑑𝑡𝜑2(𝑡)+0.5sin(𝜑4(𝑡))𝑑𝑑𝑡𝜑4(𝑡))−54.51345(2.0cos(𝜑0(𝑡))𝑑𝑑𝑡𝜑0(𝑡)+1.0cos(𝜑1(𝑡))𝑑𝑑𝑡𝜑1(𝑡))𝜑0(𝑡)𝑑𝑑𝑡𝜑0(𝑡)+54.51345(𝜑0(𝑡)+sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡)))((0.5𝑑𝑑𝑡𝜑0(𝑡)+0.5𝑑𝑑𝑡𝜑1(𝑡)−0.5𝑑𝑑𝑡𝜑2(𝑡))(𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))cos(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+(0.5𝑑2𝑑𝑡2𝜑0(𝑡)+0.5𝑑2𝑑𝑡2𝜑1(𝑡)−0.5𝑑2𝑑𝑡2𝜑2(𝑡))sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+0.5𝜑0(𝑡)𝑑2𝑑𝑡2𝜑0(𝑡)+0.5sin(𝜑1(𝑡))𝑑2𝑑𝑡2𝜑1(𝑡)+1.0sin(𝜑2(𝑡))𝑑2𝑑𝑡2𝜑2(𝑡)+0.5cos(𝜑0(𝑡))(𝑑𝑑𝑡𝜑0(𝑡))2+0.5cos(𝜑1(𝑡))(𝑑𝑑𝑡𝜑1(𝑡))2+1.0cos(𝜑2(𝑡))(𝑑𝑑𝑡𝜑2(𝑡))2)+54.51345(𝜑0(𝑡)+2.0sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡)))((𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))(1.0𝑑𝑑𝑡𝜑0(𝑡)+1.0𝑑𝑑𝑡𝜑1(𝑡)−1.0𝑑𝑑𝑡𝜑2(𝑡))cos(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+(1.0𝑑2𝑑𝑡2𝜑0(𝑡)+1.0𝑑2𝑑𝑡2𝜑1(𝑡)−1.0𝑑2𝑑𝑡2𝜑2(𝑡))sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+0.5𝜑0(𝑡)𝑑2𝑑𝑡2𝜑0(𝑡)+0.5sin(𝜑1(𝑡))𝑑2𝑑𝑡2𝜑1(𝑡)+1.0sin(𝜑2(𝑡))𝑑2𝑑𝑡2𝜑2(𝑡)+0.5sin(𝜑4(𝑡))𝑑2𝑑𝑡2𝜑4(𝑡)+0.5cos(𝜑0(𝑡))(𝑑𝑑𝑡𝜑0(𝑡))2+0.5cos(𝜑1(𝑡))(𝑑𝑑𝑡𝜑1(𝑡))2+1.0cos(𝜑2(𝑡))(𝑑𝑑𝑡𝜑2(𝑡))2+0.5cos(𝜑4(𝑡))(𝑑𝑑𝑡𝜑4(𝑡))2)+54.51345(cos(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))−2.0cos(𝜑0(𝑡)))(−(0.5𝑑𝑑𝑡𝜑0(𝑡)+0.5𝑑𝑑𝑡𝜑1(𝑡)−0.5𝑑𝑑𝑡𝜑2(𝑡))(𝑑𝑑𝑡𝜑0(𝑡)+𝑑𝑑𝑡𝜑1(𝑡)−𝑑𝑑𝑡𝜑2(𝑡))sin(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+(0.5𝑑2𝑑𝑡2𝜑0(𝑡)+0.5𝑑2𝑑𝑡2𝜑1(𝑡)−0.5𝑑2𝑑𝑡2𝜑2(𝑡))cos(𝜑0(𝑡)+𝜑1(𝑡)−𝜑2(𝑡))+1.0𝜑0(𝑡)(𝑑𝑑𝑡𝜑0(𝑡))2+1.0sin(𝜑1(𝑡))(𝑑𝑑𝑡𝜑1(𝑡))2−1.0sin(𝜑2(𝑡))(𝑑𝑑𝑡𝜑2(𝑡))2−1.0cos(𝜑0(𝑡))𝑑2𝑑𝑡2𝜑0(𝑡)−1.0cos(𝜑1(𝑡))𝑑2𝑑𝑡2𝜑1(𝑡)+1.0cos(𝜑2(𝑡))𝑑2𝑑𝑡2𝜑2(𝑡))+54.51345(0.5𝜑0(𝑡)𝑑𝑑𝑡𝜑0(𝑡)+0.5sin(𝜑1(𝑡))𝑑𝑑𝑡𝜑1(𝑡)+0.5sin(𝜑2(𝑡))𝑑𝑑𝑡𝜑2(𝑡))cos(𝜑0(𝑡))𝑑𝑑𝑡𝜑0(𝑡)−54.51345(2.0cos(𝜑0(𝑡))𝑑𝑑𝑡𝜑0(𝑡)+2.0cos(𝜑1(𝑡))𝑑𝑑𝑡𝜑1(𝑡)−1.0cos(𝜑2(𝑡))𝑑𝑑𝑡𝜑2(𝑡))𝜑0(𝑡)𝑑𝑑𝑡𝜑0(𝑡)+54.51345(−2.0𝜑0(𝑡)(𝑑𝑑𝑡𝜑0(𝑡))2−1.0sin(𝜑1(𝑡))(𝑑𝑑𝑡𝜑1(𝑡))2+2.0cos(𝜑0(𝑡))𝑑2𝑑𝑡2𝜑0(𝑡)+1.0cos(𝜑1(𝑡))𝑑2𝑑𝑡2𝜑1(𝑡))cos(𝜑0(𝑡))+54.51345(−2.0𝜑0(𝑡)(𝑑𝑑𝑡𝜑0(𝑡))2−2.0sin(𝜑1(𝑡))(𝑑𝑑𝑡𝜑1(𝑡))2+1.0sin(𝜑2(𝑡))(𝑑𝑑𝑡𝜑2(𝑡))2+2.0cos(𝜑0(𝑡))𝑑2𝑑𝑡2𝜑0(𝑡)+2.0cos(𝜑1(𝑡))𝑑2𝑑𝑡2𝜑1(𝑡)−1.0cos(𝜑2(𝑡))𝑑2𝑑𝑡2𝜑2(𝑡))cos(𝜑0(𝑡))+54.51345(0.5𝜑0(𝑡)𝑑2𝑑𝑡2𝜑0(𝑡)+0.5sin(𝜑1(𝑡))𝑑2𝑑𝑡2𝜑1(𝑡)+0.5sin(𝜑2(𝑡))𝑑2𝑑𝑡2𝜑2(𝑡)+0.5cos(𝜑0(𝑡))(𝑑𝑑𝑡𝜑0(𝑡))2+0.5cos(𝜑1(𝑡))(𝑑𝑑𝑡𝜑1(𝑡))2+0.5cos(𝜑2(𝑡))(𝑑𝑑𝑡𝜑2(𝑡))2)𝜑0(𝑡)−2123.406𝑑𝑑𝑡𝜑0(𝑡)+45.427875𝑑2𝑑𝑡2𝜑0(𝑡)−2123.406𝑑𝑑𝑡𝜑1(𝑡)+9.085575𝑑2𝑑𝑡2𝜑1(𝑡)−9.085575𝑑2𝑑𝑡2𝜑2(𝑡)

和lambdify输出：

NameError                                 Traceback (most recent call last)
<ipython-input-14-ee077b324a2e> in <module>
      2 en2 = sym.lambdify([kot_0,kot_1,kot_2],f1)
      3 en3 = sym.lambdify([kot_0,kot_1,kot_2],f2)
----> 4 en1(kot_0,kot_1,kot_2)
      5 
      6 

<lambdifygenerated-4> in _lambdifygenerated(_Dummy_227, _Dummy_226, _Dummy_225)
      9   # Derivative
     10   # Derivative
---> 11 22.722525*_Dummy_227**2*Derivative(_Dummy_227, (t, 2)) + 90.8901*_Dummy_227*(-0.5*cos(_Dummy_226)*Derivative(_Dummy_226, t) - 1.0*cos(_Dummy_227)*Derivative(_Dummy_227, t))*Derivative(_Dummy_227, t) + 90.8901*_Dummy_227*(0.5*cos(_Dummy_225)*Derivative(_Dummy_225, t) - 1.0*cos(_Dummy_226)*Derivative(_Dummy_226, t) - 1.0*cos(_Dummy_227)*Derivative(_Dummy_227, t))*Derivative(_Dummy_227, t) - 22.722525*_Dummy_227*(-_Dummy_227*Derivative(_Dummy_227, (t, 2)) - sin(_Dummy_225)*Derivative(_Dummy_225, (t, 2)) - sin(_Dummy_226)*Derivative(_Dummy_226, (t, 2)) - cos(_Dummy_225)*Derivative(_Dummy_225, t)**2 - cos(_Dummy_226)*Derivative(_Dummy_226, t)**2 - cos(_Dummy_227)*Derivative(_Dummy_227, t)**2) + 0.5*Jm*(-2*Derivative(_Dummy_225, (t, 2)) + 2*Derivative(_Dummy_226, (t, 2)) + 2*Derivative(_Dummy_227, (t, 2))) + 1.0*Jm*Derivative(_Dummy_227, (t, 2)) + 45.44505*(-1.0*_Dummy_227 - 2.0*sin(-_Dummy_225 + _Dummy_226 + _Dummy_227))*(-0.5*_Dummy_227*Derivative(_Dummy_227, (t, 2)) + (-Derivative(_Dummy_225, t) + Derivative(_Dummy_226, t) + Derivative(_Dummy_227, t))*(1.0*Derivative(_Dummy_225, t) - 1.0*Derivative(_Dummy_226, t) - 1.0*Derivative(_Dummy_227, t))*cos(-_Dummy_225 + _Dummy_226 + _Dummy_227) + (1.0*Derivative(_Dummy_225, (t, 2)) - 1.0*Derivative(_Dummy_226, (t, 2)) - 1.0*Derivative(_Dummy_227, (t, 2)))*sin(-_Dummy_225 + _Dummy_226 + _Dummy_227) - 1.0*sin(_Dummy_225)*Derivative(_Dummy_225, (t, 2)) - 0.5*sin(_Dummy_226)*Derivative(_Dummy_226, (t, 2)) - 0.5*sin(varphi_4(t))*Der

https://imgur.com/a/2UOW0NR

EDIT2：经过一段时间的简化，我得到了3个常微分方程式。但是它们是渐近形式的，我怎么能用数字来求解呢？

−800000.0𝜑0+800000.0𝜑2−1770.174𝜑˙0+242.3736𝜑¨0−1770.174𝜑˙1+166.63185𝜑¨1−75.74175𝜑¨2+4245.8661 

−1200000.0𝜑1+400000.0𝜑2−1770.174𝜑˙0+166.63185𝜑¨0−1770.174𝜑˙1+151.4835𝜑¨1−75.74175𝜑¨2+2830.5774 

800000.0𝜑0+400000.0𝜑1−2000000.0𝜑2−75.74175𝜑¨0−75.74175𝜑¨1+60.5934𝜑¨2−1415.2887

推荐答案

我将以二级摆为例介绍欧拉-拉格朗日公式的执行，这是一个不平凡的、完整的、标准的例子，这不需要使用渐近物理中的专门函数，只使用基本的渐近微分和代码编写工具。希望它足够通用，因为第二部分应该是与问题无关的，因此也直接适用于您的情况。

from sympy import sin, cos, Symbol, symbols, solve
from sympy.utilities.lambdify import lambdify

物理模型的拉格朗日

首先，使用两个角作为主要依赖函数，建立拉格朗日的物理设置，并通过笛卡尔坐标构造动能和势能。

# Symbols for the parameters of the problem
t,m1,m2,l1,l2,g = symbols("t,m_1 m_2 l_1 l_2 g")
# Variables of the problem
th1, th2 = Function("θ_1")(t), Function("θ_2")(t)
x1,y1 = l1*sin(th1), -l1*cos(th1)
x2,y2 = x1+l2*sin(th2), y1-l2*cos(th2)

# kinetic energy
vx1,vy1,vx2,vy2 = ( xx.diff(t) for xx in (x1,y1,x2,y2))
K1 = m1/2*(vx1**2+vy1**2)
K2 = m2/2*(vx2**2+vy2**2)
K = K1+K2
# potential energy
V = g*(m1*y1+m2*y2)
# Lagrangian
L = K - V
L = L.expand().simplify()

使用抽象参数数组和坐标向量获取抽象处理

params = [m1, l1, m2, l2, g]
q = [th1, th2]
dotq = [ qq.diff(t) for qq in q]

从欧拉-拉格朗日到一阶常微分方程组

如果到角度的二阶导数，结果表达式会变得相当混乱，得到更结构化的方法，并希望使用脉冲变量进行相应的更快计算

pk = diff(L, dotqk)
d(pk)/dt = diff(L, qk)

其中第一个关系被视为要从p计算dotq的方程系统。

准备新变量，准备用简单变量替换函数符号。

N = len(q)
p = [ Symbol(f"p_{k+1}") for k in range(N)]
dotq_func, dotq = dotq,  [ Symbol(f"Dq_{k+1}") for k in range(N)]
q_func, q = q, [ Symbol(f"q_{k+1}") for k in range(N)]

现在将所有函数项替换为简单变量

L=L.subs(list(zip(dotq_func, dotq))).subs(list(zip(q_func, q)))

现在建立要从q,p

计算dotq的函数

p_eqns = [ p[k] - L.diff(dotq[k]) for k in range(N)]
dotq_expr = solve(p_eqns, dotq)
dotq_func = lambdify([*q, *p, *params],[ dotq_expr[dq] for dq in dotq])

接下来生成一个函数，用于计算p

的导数

dotp_expr = [ L.diff(q[k]) for k in range(N)]
dotp_func = lambdify([*q,*dotq,*params],dotp_expr)

现在将生成的函数组装成一个完整的ODE函数，并用它解决一个测试问题以确认它可以工作

def odefunc(t,u,args):
    q, p = u[:N], u[N:]
    dotq = dotq_func(*q, *p, *args)
    dotp = dotp_func(*q, *dotq, *args)
    return [*dotq, *dotp]

通过数值解确认

myparams = [1,10,5,5,9.81]
t = np.linspace(0,25,301)

u = odeint(odefunc,[2,1.2, 0,0],t,args=(myparams,), tfirst=True)

%matplotlib inline
plt.figure(figsize=(8,5))
plt.plot(t,u[:,0],t,u[:,1]); 
plt.grid(); plt.show()

这将为角度函数生成以下曲线图

Outlook

如果动力学项(如经典力学中常见的那样)保证速度是二次项，则可能会有更好的性能。然后，通过直接提取这种二次型的矩阵，可以将脉冲到速度转换中的系统解委托给数值线性系统求解器，从而不保留用于矩阵求逆的符号表达式。这些在更高维度中可能很大。

这篇关于从渐近性出发求解二阶微分方程组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！