本文介绍了从渐近性出发求解二阶微分方程组的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我正在做一个多自由度动力学问题，使用的是二阶拉格朗日方程。我用渐近解出了运动方程。现在这些方程在计算导数后变得相当长，尽管符号简化不能进一步简化它。我的问题实际上是如何从这里解决这个三个二阶颂歌的系统。我不知道如何转换这些公式，这样它们才能与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
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
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
如果动力学项(如经典力学中常见的那样)保证速度是二次项，则可能会有更好的性能。然后，通过直接提取这种二次型的矩阵，可以将脉冲到速度转换中的系统解委托给数值线性系统求解器，从而不保留用于矩阵求逆的符号表达式。这些在更高维度中可能很大。
这篇关于从渐近性出发求解二阶微分方程组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！