在 sympy 中解决的性能问题

Question

第一次发帖。我目前正在尝试在 Python 中使用 sympy 为 17 个变量求解 17 个符号方程。我的方程在 17 个变量中是线性的。我遇到的问题是求解函数非常慢 - 我已经让我的程序运行了两个多小时，但它仍然没有产生结果。

我尝试将标志“简化”设置为 false，将“理性”设置为 false，但我的程序仍然需要很长时间才能运行（几个小时后它没有产生结果）。我也尝试过使用linsolve，但我让它运行一夜后我的程序没有完成。

我还将我的方程插入到 MATLAB 中，将它们转换为矩阵形式，并使用 MATLAB 的反斜杠命令求解它们。 MATLAB 在 2 分钟多一点的时间内给出了我的结果。

有人对如何提高我的 Python 代码的性能有任何建议吗？我希望能够在 Python 中完成所有工作，而不必依赖 MATLAB 的符号工具箱。

Python 代码：

import sympy as sym
from sympy import cos as cos
from sympy import sin as sin
import time

#Starting timer
t0 = time.time()

#Defining symbolic variables
x1,y1,theta1,x1dot,y1dot,theta1dot,x1ddot,y1ddot,theta1ddot,x2,y2,theta2,x2dot,y2dot,theta2dot, \
    x2ddot,y2ddot,theta2ddot,x3,y3,theta3,x3dot,y3dot,theta3dot,x3ddot,y3ddot,theta3ddot, \
    R1x,R1y,J1x,J1y,J2x,J2y,R2x,R2y,l1,l2,l3,m1,m2,m3,g = sym.symbols("x1 y1 theta1 x1dot y1dot theta1dot x1ddot y1ddot "
    "theta1ddot x2 y2 theta2 x2dot y2dot theta2dot x2ddot y2ddot theta2ddot x3 y3 theta3 x3dot y3dot theta3dot x3ddot "
    "y3ddot theta3ddot R1x R1y J1x J1y J2x J2y R2x R2y l1 l2 l3 m1 m2 m3 g",real=True)

#Defining parameters
d1 = l1/2; d2 = l2/2; d3 = l3/2
I1 = (1/12)*m1*l1**2
I2 = (1/12)*m2*l2**2
I3 = (1/12)*m3*l3**2

#Defining unit vectors
ihat = sym.Matrix([1,0,0])
jhat = sym.Matrix([0,1,0])
khat = sym.Matrix([0,0,1])

#Defining useful position vectors
rg1 = d1*sin(theta1)*ihat - d1*cos(theta1)*jhat
rp1 = l1*sin(theta1)*ihat - l1*cos(theta1)*jhat
rg2p1 = d2*sin(theta2)*ihat - d2*cos(theta2)*jhat
rp2p1 = l2*sin(theta2)*ihat - l2*cos(theta2)*jhat
rg3p2 = d3*sin(theta3)*ihat - d3*cos(theta3)*jhat
rp3p2 = l3*sin(theta3)*ihat - l3*cos(theta3)*jhat

#Defining useful acceleration vectors
ap1 = (theta1ddot*khat).cross(rp1) + (theta1dot*khat).cross((theta1dot*khat).cross(rp1))
ap2 = ap1 + (theta2ddot*khat).cross(rp2p1) + (theta2dot*khat).cross((theta2dot*khat).cross(rp2p1))

#Defining equations
eqn1 = sym.Eq( m1*x1ddot, R1x + J1x)
eqn2 = sym.Eq( m1*y1ddot, -m1*g + R1y - J1y)
eqn3 = sym.Eq( m2*x2ddot, -J1x + J2x)
eqn4 = sym.Eq( m2*y2ddot, -m2*g + J1y - J2y)
eqn5 = sym.Eq( m3*x3ddot, -J2x + R2x)
eqn6 = sym.Eq( m3*y3ddot, -m3*g + J2y + R2y)
eqn7 = sym.Eq( (rg1.cross(-m1*g*jhat) + rp1.cross(J1x*ihat - J1y*jhat))[2],
               (rg1.cross(m1*(x1ddot*ihat + y1ddot*jhat)) + I1*theta1ddot*khat)[2])
eqn8 = sym.Eq( (rg2p1.cross(-m2*g*jhat) + rp2p1.cross(J2x*ihat - J2y*jhat))[2],
               (rg2p1.cross(m2*(x2ddot*ihat + y2ddot*jhat)) + I2*theta2ddot*khat)[2])
eqn9 = sym.Eq( (rg3p2.cross(-m3*g*jhat) + rp3p2.cross(R2x*ihat + R2y*jhat))[2],
               (rg3p2.cross(m3*(x3ddot*ihat + y3ddot*jhat)) + I3*theta3ddot*khat)[2])
eqn10 = sym.Eq( (x1ddot*ihat + y1ddot*jhat)[0],
                ((theta1ddot*khat).cross(rg1) + (theta1dot*khat).cross((theta1dot*khat).cross(rg1)))[0])
eqn11 = sym.Eq( (x1ddot*ihat + y1ddot*jhat)[1],
                ((theta1ddot*khat).cross(rg1) + (theta1dot*khat).cross((theta1dot*khat).cross(rg1)))[1])
eqn12 = sym.Eq( (x2ddot*ihat + y2ddot*jhat)[0],
                (ap1 + (theta2ddot*khat).cross(rg2p1) + (theta2dot*khat).cross((theta2dot*khat).cross(rg2p1)))[0])
eqn13 = sym.Eq( (x2ddot*ihat + y2ddot*jhat)[1],
                (ap1 + (theta2ddot*khat).cross(rg2p1) + (theta2dot*khat).cross((theta2dot*khat).cross(rg2p1)))[1])
eqn14 = sym.Eq( (x3ddot*ihat + y3ddot*jhat)[0],
                (ap2 + (theta3ddot*khat).cross(rg3p2) + (theta3dot*khat).cross((theta3dot*khat).cross(rg3p2)))[0])
eqn15 = sym.Eq( (x3ddot*ihat + y3ddot*jhat)[1],
                (ap2 + (theta3ddot*khat).cross(rg3p2) + (theta3dot*khat).cross((theta3dot*khat).cross(rg3p2)))[1])
eqn16 = sym.Eq( (ap2 + (theta3ddot*khat).cross(rp3p2) + (theta3dot*khat).cross((theta3dot*khat).cross(rp3p2)))[0],0)
eqn17 = sym.Eq( (ap2 + (theta3ddot*khat).cross(rp3p2) + (theta3dot*khat).cross((theta3dot*khat).cross(rp3p2)))[1],0)

#Lists of equations and variables
eqns = [eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,eqn11,eqn12,eqn13,eqn14,eqn15,eqn16,eqn17]
vars = [x1ddot,y1ddot,theta1ddot,x2ddot,y2ddot,theta2ddot,x3ddot,y3ddot,theta3ddot,R1x,R1y,J1x,J1y,J2x,J2y,R2x,R2y]

#Generating u vector
u = sym.solve(eqns,vars,simplify=False,rational=False)

elapsedTime = time.time()-t0

print(u)
print(type(u))
print(elapsedTime, " seconds")

MATLAB 代码：

clear; close all; clc;
tic
%Defining variables
syms x1 y1 theta1 x1dot y1dot theta1dot x1ddot y1ddot theta1ddot x2 y2 theta2 x2dot y2dot theta2dot ...
    x2ddot y2ddot theta2ddot x3 y3 theta3 x3dot y3dot theta3dot x3ddot y3ddot theta3ddot ...
    l1 l2 l3 g R1x R1y J1x J1y J2x J2y R2x R2y m1 m2 m3 real

d1 = l1/2; d2 = l2/2; d3 = l3/2;

%Defining moments of inertia
I1 = (1/12)*m1*l1^2;
I2 = (1/12)*m2*l2^2;
I3 = (1/12)*m3*l3^2;

%Defining unit vectors
ihat = [1 0 0];
jhat = [0 1 0];
khat = [0 0 1];

%Defining useful position vectors
rg1 = d1*sin(theta1)*ihat - d1*cos(theta1)*jhat;
rp1 = l1*sin(theta1)*ihat - l1*cos(theta1)*jhat;
rg2p1 = d2*sin(theta2)*ihat - d2*cos(theta2)*jhat;
rp2p1 = l2*sin(theta2)*ihat - l2*cos(theta2)*jhat;
rg3p2 = d3*sin(theta3)*ihat - d3*cos(theta3)*jhat;
rp3p2 = l3*sin(theta3)*ihat - l3*cos(theta3)*jhat;

%Defining useful acceleration vectors
ap1 = cross(theta1ddot*khat,rp1) + cross(theta1dot*khat,cross(theta1dot*khat,rp1));
ap2 = ap1 + cross(theta2ddot*khat,rp2p1) + cross(theta2dot*khat,cross(theta2dot*khat,rp2p1));

%Defining equations of motion
eq1 = m1*x1ddot == R1x + J1x;
eq2 = m1*y1ddot == -m1*g + R1y - J1y;
eq3 = m2*x2ddot == -J1x + J2x;
eq4 = m2*y2ddot == -m2*g + J1y - J2y;
eq5 = m3*x3ddot == -J2x + R2x;
eq6 = m3*y3ddot == -m3*g + J2y + R2y;
eq7 = (cross(rg1,-m1*g*jhat) + cross(rp1,J1x*ihat - J1y*jhat)) == ...
    (cross(rg1,m1*(x1ddot*ihat + y1ddot*jhat)) + I1*theta1ddot*khat);
eq7 = eq7(3);
eq8 = cross(rg2p1,-m2*g*jhat) + cross(rp2p1,J2x*ihat-J2y*jhat) == ...
    cross(rg2p1,m2*(x2ddot*ihat+y2ddot*jhat)) + I2*theta2ddot*khat;
eq8 = eq8(3);
eq9 = cross(rg3p2,-m3*g*jhat) + cross(rp3p2,R2x*ihat+R2y*jhat) == ...
    cross(rg3p2,m3*(x3ddot*ihat+y3ddot*jhat))+I3*theta3ddot*khat;
eq9 = eq9(3);
eqns10_11 = x1ddot*ihat + y1ddot*jhat == ...
    cross(theta1ddot*khat,rg1) + cross(theta1dot*khat,cross(theta1dot*khat,rg1));
eq10 = eqns10_11(1);
eq11 = eqns10_11(2);
eqns12_13 = x2ddot*ihat + y2ddot*jhat == ...
    ap1 + cross(theta2ddot*khat,rg2p1) + cross(theta2dot*khat,cross(theta2dot*khat,rg2p1));
eq12 = eqns12_13(1);
eq13 = eqns12_13(2);
eqns14_15 = x3ddot*ihat + y3ddot*jhat == ...
    ap2 + cross(theta3ddot*khat,rg3p2) + cross(theta3dot*khat,cross(theta3dot*khat,rg3p2));
eq14 = eqns14_15(1);
eq15 = eqns14_15(2);
eqns16_17 = ap2 + cross(theta3ddot*khat,rp3p2) + cross(theta3dot*khat,cross(theta3dot*khat,rp3p2)) == 0;
eq16 = eqns16_17(1);
eq17 = eqns16_17(2);

%Putting equations in matrix form
eqns = [eq1 eq2 eq3 eq4 eq5 eq6 eq7 eq8 eq9 eq10 eq11 eq12 eq13 eq14 eq15 eq16 eq17];
vars = [x1ddot y1ddot theta1ddot x2ddot y2ddot theta2ddot x3ddot y3ddot theta3ddot ...
    R1x R1y J1x J1y J2x J2y R2x R2y];
[A,b] = equationsToMatrix(eqns,vars);

u = A\b;

toc

MATLAB 输出（这只是输出的一部分。我的帖子正文限制为 30000 个字符，而整个输出发布后，帖子包含几乎 67000 个字符。我做了一些检查，这个输出是正确。它也很大——它是一个 17 x 1 的向量，每个条目用换行符分隔。我已经粘贴了前几个条目）

(3*g*m2*sin(2*theta2) - 9*g*m2*sin(2*theta1) - 6*g*m1*sin(2*theta1) - 3*g*m3*sin(2*theta1) - 3*g*m2*sin(2*theta3) + 3*g*m3*sin(2*theta2) - 3*g*m3*sin(2*theta3) + 3*g*m1*sin(2*theta1 - 2*theta2 + 2*theta3) + 3*g*m1*sin(2*theta1 + 2*theta2 - 2*theta3) + 6*g*m2*sin(2*theta1 - 2*theta2 + 2*theta3) + 3*g*m2*sin(2*theta1 + 2*theta2 - 2*theta3) + 3*g*m3*sin(2*theta1 - 2*theta2 + 2*theta3) - 3*g*m2*sin(2*theta1 - 2*theta2) + 3*g*m2*sin(2*theta1 - 2*theta3) - 3*g*m3*sin(2*theta1 - 2*theta2) - 3*g*m2*sin(2*theta2 - 2*theta3) + 3*g*m3*sin(2*theta1 - 2*theta3) - 3*g*m3*sin(2*theta2 - 2*theta3) - 12*l1*m2*theta1dot^2*sin(theta1 - 2*theta2) + 4*l1*m2*theta1dot^2*sin(theta1 - 2*theta3) - 8*l1*m3*theta1dot^2*sin(theta1 - 2*theta2) + 2*l2*m2*theta2dot^2*sin(theta2 - 2*theta3) - 10*l2*m2*theta2dot^2*sin(2*theta1 - theta2) - 8*l2*m3*theta2dot^2*sin(2*theta1 - theta2) - 2*l3*m2*theta3dot^2*sin(2*theta1 - theta3) + 6*l3*m2*theta3dot^2*sin(2*theta2 - theta3) - 4*l3*m3*theta3dot^2*sin(2*theta1 - theta3) + 4*l3*m3*theta3dot^2*sin(2*theta2 - theta3) - 8*l1*m1*theta1dot^2*sin(theta1) - 20*l1*m2*theta1dot^2*sin(theta1) - 8*l1*m3*theta1dot^2*sin(theta1) + 10*l2*m2*theta2dot^2*sin(theta2) + 8*l2*m3*theta2dot^2*sin(theta2) + 2*l3*m2*theta3dot^2*sin(theta3) + 4*l3*m3*theta3dot^2*sin(theta3) + 4*l1*m1*theta1dot^2*sin(theta1 - 2*theta2 + 2*theta3) + 4*l1*m1*theta1dot^2*sin(theta1 + 2*theta2 - 2*theta3) + 6*l1*m2*theta1dot^2*sin(theta1 - 2*theta2 + 2*theta3) + 6*l1*m2*theta1dot^2*sin(theta1 + 2*theta2 - 2*theta3) + 2*l2*m2*theta2dot^2*sin(2*theta1 + theta2 - 2*theta3) - 6*l3*m2*theta3dot^2*sin(2*theta1 - 2*theta2 + theta3) - 4*l3*m3*theta3dot^2*sin(2*theta1 - 2*theta2 + theta3))/(8*(2*m1 + 5*m2 + 2*m3 - 3*m2*cos(2*theta1 - 2*theta2) - 2*m1*cos(2*theta2 - 2*theta3) + m2*cos(2*theta1 - 2*theta3) - 2*m3*cos(2*theta1 - 2*theta2) - 3*m2*cos(2*theta2 - 2*theta3)))

-(6*g*m1 + 9*g*m2 + 3*g*m3 - 6*g*m1*cos(2*theta1) - 9*g*m2*cos(2*theta1) + 3*g*m2*cos(2*theta2) - 3*g*m3*cos(2*theta1) - 3*g*m2*cos(2*theta3) + 3*g*m3*cos(2*theta2) - 3*g*m3*cos(2*theta3) + 3*g*m1*cos(2*theta1 - 2*theta2 + 2*theta3) + 3*g*m1*cos(2*theta1 + 2*theta2 - 2*theta3) + 6*g*m2*cos(2*theta1 - 2*theta2 + 2*theta3) + 3*g*m2*cos(2*theta1 + 2*theta2 - 2*theta3) + 3*g*m3*cos(2*theta1 - 2*theta2 + 2*theta3) - 3*g*m2*cos(2*theta1 - 2*theta2) - 6*g*m1*cos(2*theta2 - 2*theta3) + 3*g*m2*cos(2*theta1 - 2*theta3) - 3*g*m3*cos(2*theta1 - 2*theta2) - 9*g*m2*cos(2*theta2 - 2*theta3) + 3*g*m3*cos(2*theta1 - 2*theta3) - 3*g*m3*cos(2*theta2 - 2*theta3) + 12*l1*m2*theta1dot^2*cos(theta1 - 2*theta2) - 4*l1*m2*theta1dot^2*cos(theta1 - 2*theta3) + 8*l1*m3*theta1dot^2*cos(theta1 - 2*theta2) - 2*l2*m2*theta2dot^2*cos(theta2 - 2*theta3) - 10*l2*m2*theta2dot^2*cos(2*theta1 - theta2) - 8*l2*m3*theta2dot^2*cos(2*theta1 - theta2) - 2*l3*m2*theta3dot^2*cos(2*theta1 - theta3) + 6*l3*m2*theta3dot^2*cos(2*theta2 - theta3) - 4*l3*m3*theta3dot^2*cos(2*theta1 - theta3) + 4*l3*m3*theta3dot^2*cos(2*theta2 - theta3) - 8*l1*m1*theta1dot^2*cos(theta1) - 20*l1*m2*theta1dot^2*cos(theta1) - 8*l1*m3*theta1dot^2*cos(theta1) + 10*l2*m2*theta2dot^2*cos(theta2) + 8*l2*m3*theta2dot^2*cos(theta2) + 2*l3*m2*theta3dot^2*cos(theta3) + 4*l3*m3*theta3dot^2*cos(theta3) + 4*l1*m1*theta1dot^2*cos(theta1 - 2*theta2 + 2*theta3) + 4*l1*m1*theta1dot^2*cos(theta1 + 2*theta2 - 2*theta3) + 6*l1*m2*theta1dot^2*cos(theta1 - 2*theta2 + 2*theta3) + 6*l1*m2*theta1dot^2*cos(theta1 + 2*theta2 - 2*theta3) + 2*l2*m2*theta2dot^2*cos(2*theta1 + theta2 - 2*theta3) - 6*l3*m2*theta3dot^2*cos(2*theta1 - 2*theta2 + theta3) - 4*l3*m3*theta3dot^2*cos(2*theta1 - 2*theta2 + theta3))/(8*(2*m1 + 5*m2 + 2*m3 - 3*m2*cos(2*theta1 - 2*theta2) - 2*m1*cos(2*theta2 - 2*theta3) + m2*cos(2*theta1 - 2*theta3) - 2*m3*cos(2*theta1 - 2*theta2) - 3*m2*cos(2*theta2 - 2*theta3)))

-(6*g*m1*sin(theta1) + 9*g*m2*sin(theta1) + 3*g*m3*sin(theta1) - 3*g*m1*sin(theta1 - 2*theta2 + 2*theta3) - 3*g*m1*sin(theta1 + 2*theta2 - 2*theta3) - 6*g*m2*sin(theta1 - 2*theta2 + 2*theta3) - 3*g*m2*sin(theta1 + 2*theta2 - 2*theta3) - 3*g*m3*sin(theta1 - 2*theta2 + 2*theta3) + 3*g*m2*sin(theta1 - 2*theta2) - 3*g*m2*sin(theta1 - 2*theta3) + 3*g*m3*sin(theta1 - 2*theta2) - 3*g*m3*sin(theta1 - 2*theta3) + 10*l2*m2*theta2dot^2*sin(theta1 - theta2) + 8*l2*m3*theta2dot^2*sin(theta1 - theta2) + 2*l3*m2*theta3dot^2*sin(theta1 - theta3) + 4*l3*m3*theta3dot^2*sin(theta1 - theta3) + 6*l1*m2*theta1dot^2*sin(2*theta1 - 2*theta2) - 2*l1*m2*theta1dot^2*sin(2*theta1 - 2*theta3) + 4*l1*m3*theta1dot^2*sin(2*theta1 - 2*theta2) - 2*l2*m2*theta2dot^2*sin(theta1 + theta2 - 2*theta3) + 6*l3*m2*theta3dot^2*sin(theta1 - 2*theta2 + theta3) + 4*l3*m3*theta3dot^2*sin(theta1 - 2*theta2 + theta3))/(2*l1*(2*m1 + 5*m2 + 2*m3 - 3*m2*cos(2*theta1 - 2*theta2) - 2*m1*cos(2*theta2 - 2*theta3) + m2*cos(2*theta1 - 2*theta3) - 2*m3*cos(2*theta1 - 2*theta2) - 3*m2*cos(2*theta2 - 2*theta3)))

-(9*g*m1*sin(2*theta1) - 3*g*m1*sin(2*theta2) + 12*g*m2*sin(2*theta1) + 3*g*m1*sin(2*theta3) - 6*g*m2*sin(2*theta2) + 3*g*m3*sin(2*theta1) + 12*g*m2*sin(2*theta3) - 3*g*m3*sin(2*theta2) + 9*g*m3*sin(2*theta3) - 6*g*m1*sin(2*theta1 - 2*theta2 + 2*theta3) - 3*g*m1*sin(2*theta1 + 2*theta2 - 2*theta3) - 3*g*m2*sin(2*theta2 - 2*theta1 + 2*theta3) - 12*g*m2*sin(2*theta1 - 2*theta2 + 2*theta3) - 3*g*m2*sin(2*theta1 + 2*theta2 - 2*theta3) - 3*g*m3*sin(2*theta2 - 2*theta1 + 2*theta3) - 6*g*m3*sin(2*theta1 - 2*theta2 + 2*theta3) - 3*g*m1*sin(2*theta1 - 2*theta2) + 3*g*m1*sin(2*theta1 - 2*theta3) - 3*g*m1*sin(2*theta2 - 2*theta3) + 3*g*m3*sin(2*theta1 - 2*theta2) - 3*g*m3*sin(2*theta1 - 2*theta3) + 3*g*m3*sin(2*theta2 - 2*theta3) - 4*l1*m1*theta1dot^2*sin(theta1 - 2*theta2) + 4*l1*m1*theta1dot^2*sin(theta1 - 2*theta3) + 6*l1*m2*theta1dot^2*sin(theta1 - 2*theta2) - 2*l1*m2*theta1dot^2*sin(theta1 - 2*theta3) + 8*l1*m3*theta1dot^2*sin(theta1 - 2*theta2) + 8*l2*m1*theta2dot^2*sin(theta2 - 2*theta3) + 8*l2*m2*theta2dot^2*sin(theta2 - 2*theta3) + 8*l2*m2*theta2dot^2*sin(2*theta1 - theta2) + 8*l2*m3*theta2dot^2*sin(2*theta1 - theta2) + 8*l3*m1*theta3dot^2*sin(2*theta2 - theta3) - 2*l3*m2*theta3dot^2*sin(2*theta1 - theta3) + 6*l3*m2*theta3dot^2*sin(2*theta2 - theta3) + 4*l3*m3*theta3dot^2*sin(2*theta1 - theta3) - 4*l3*m3*theta3dot^2*sin(2*theta2 - theta3) + 12*l1*m1*theta1dot^2*sin(theta1) + 22*l1*m2*theta1dot^2*sin(theta1) + 8*l1*m3*theta1dot^2*sin(theta1) + 8*l2*m1*theta2dot^2*sin(theta2) - 8*l2*m3*theta2dot^2*sin(theta2) - 8*l3*m1*theta3dot^2*sin(theta3) - 22*l3*m2*theta3dot^2*sin(theta3) - 12*l3*m3*theta3dot^2*sin(theta3) - 8*l1*m1*theta1dot^2*sin(theta1 - 2*theta2 + 2*theta3) - 4*l1*m1*theta1dot^2*sin(theta1 + 2*theta2 - 2*theta3) - 12*l1*m2*theta1dot^2*sin(theta1 - 2*theta2 + 2*theta3) - 6*l1*m2*theta1dot^2*sin(theta1 + 2*theta2 - 2*theta3) + 2*l2*m2*theta2dot^2*sin(theta2 - 2*theta1 + 2*theta3) - 2*l2*m2*theta2dot^2*sin(2*theta1 + theta2 - 2*theta3) + 6*l3*m2*theta3dot^2*sin(2*theta2 - 2*theta1 + theta3) + 12*l3*m2*theta3dot^2*sin(2*theta1 - 2*theta2 + theta3) + 4*l3*m3*theta3dot^2*sin(2*theta2 - 2*theta1 + theta3) + 8*l3*m3*theta3dot^2*sin(2*theta1 - 2*theta2 + theta3))/(8*(2*m1 + 5*m2 + 2*m3 - 3*m2*cos(2*theta1 - 2*theta2) - 2*m1*cos(2*theta2 - 2*theta3) + m2*cos(2*theta1 - 2*theta3) - 2*m3*cos(2*theta1 - 2*theta2) - 3*m2*cos(2*theta2 - 2*theta3)))

-(9*g*m1 + 18*g*m2 + 9*g*m3 - 9*g*m1*cos(2*theta1) + 3*g*m1*cos(2*theta2) - 12*g*m2*cos(2*theta1) - 3*g*m1*cos(2*theta3) + 6*g*m2*cos(2*theta2) - 3*g*m3*cos(2*theta1) - 12*g*m2*cos(2*theta3) + 3*g*m3*cos(2*theta2) - 9*g*m3*cos(2*theta3) + 6*g*m1*cos(2*theta1 - 2*theta2 + 2*theta3) + 3*g*m1*cos(2*theta1 + 2*theta2 - 2*theta3) + 3*g*m2*cos(2*theta2 - 2*theta1 + 2*theta3) + 12*g*m2*cos(2*theta1 - 2*theta2 + 2*theta3) + 3*g*m2*cos(2*theta1 + 2*theta2 - 2*theta3) + 3*g*m3*cos(2*theta2 - 2*theta1 + 2*theta3) + 6*g*m3*cos(2*theta1 - 2*theta2 + 2*theta3) - 3*g*m1*cos(2*theta1 - 2*theta2) + 3*g*m1*cos(2*theta1 - 2*theta3) - 12*g*m2*cos(2*theta1 - 2*theta2) - 9*g*m1*cos(2*theta2 - 2*theta3) + 6*g*m2*cos(2*theta1 - 2*theta3) - 9*g*m3*cos(2*theta1 - 2*theta2) - 12*g*m2*cos(2*theta2 - 2*theta3) + 3*g*m3*cos(2*theta1 - 2*theta3) - 3*g*m3*cos(2*theta2 - 2*theta3) - 4*l1*m1*theta1dot^2*cos(theta1 - 2*theta2) + 4*l1*m1*theta1dot^2*cos(theta1 - 2*theta3) + 6*l1*m2*theta1dot^2*cos(theta1 - 2*theta2) - 2*l1*m2*theta1dot^2*cos(theta1 - 2*theta3) + 8*l1*m3*theta1dot^2*cos(theta1 - 2*theta2) + 8*l2*m1*theta2dot^2*cos(theta2 - 2*theta3) + 8*l2*m2*theta2dot^2*cos(theta2 - 2*theta3) - 8*l2*m2*theta2dot^2*cos(2*theta1 - theta2) - 8*l2*m3*theta2dot^2*cos(2*theta1 - theta2) - 8*l3*m1*theta3dot^2*cos(2*theta2 - theta3) + 2*l3*m2*theta3dot^2*cos(2*theta1 - theta3) - 6*l3*m2*theta3dot^2*cos(2*theta2 - theta3) - 4*l3*m3*theta3dot^2*cos(2*theta1 - theta3) + 4*l3*m3*theta3dot^2*cos(2*theta2 - theta3) - 12*l1*m1*theta1dot^2*cos(theta1) - 22*l1*m2*theta1dot^2*cos(theta1) - 8*l1*m3*theta1dot^2*cos(theta1) - 8*l2*m1*theta2dot^2*cos(theta2) + 8*l2*m3*theta2dot^2*cos(theta2) + 8*l3*m1*theta3dot^2*cos(theta3) + 22*l3*m2*theta3dot^2*cos(theta3) + 12*l3*m3*theta3dot^2*cos(theta3) + 8*l1*m1*theta1dot^2*cos(theta1 - 2*theta2 + 2*theta3) + 4*l1*m1*theta1dot^2*cos(theta1 + 2*theta2 - 2*theta3) + 12*l1*m2*theta1dot^2*cos(theta1 - 2*theta2 + 2*theta3) + 6*l1*m2*theta1dot^2*cos(theta1 + 2*theta2 - 2*theta3) - 2*l2*m2*theta2dot^2*cos(theta2 - 2*theta1 + 2*theta3) + 2*l2*m2*theta2dot^2*cos(2*theta1 + theta2 - 2*theta3) - 6*l3*m2*theta3dot^2*cos(2*theta2 - 2*theta1 + theta3) - 12*l3*m2*theta3dot^2*cos(2*theta1 - 2*theta2 + theta3) - 4*l3*m3*theta3dot^2*cos(2*theta2 - 2*theta1 + theta3) - 8*l3*m3*theta3dot^2*cos(2*theta1 - 2*theta2 + theta3))/(8*(2*m1 + 5*m2 + 2*m3 - 3*m2*cos(2*theta1 - 2*theta2) - 2*m1*cos(2*theta2 - 2*theta3) + m2*cos(2*theta1 - 2*theta3) - 2*m3*cos(2*theta1 - 2*theta2) - 3*m2*cos(2*theta2 - 2*theta3)))

(3*g*m1*sin(theta2) - 3*g*m3*sin(theta2) - 3*g*m1*sin(2*theta1 + theta2 - 2*theta3) + 3*g*m2*sin(theta2 - 2*theta1 + 2*theta3) - 3*g*m2*sin(2*theta1 + theta2 - 2*theta3) + 3*g*m3*sin(theta2 - 2*theta1 + 2*theta3) + 3*g*m1*sin(theta2 - 2*theta3) + 6*g*m2*sin(theta2 - 2*theta3) + 3*g*m3*sin(theta2 - 2*theta3) + 3*g*m1*sin(2*theta1 - theta2) + 6*g*m2*sin(2*theta1 - theta2) + 3*g*m3*sin(2*theta1 - theta2) + 4*l1*m1*theta1dot^2*sin(theta1 - theta2) + 18*l1*m2*theta1dot^2*sin(theta1 - theta2) + 8*l1*m3*theta1dot^2*sin(theta1 - theta2) - 8*l3*m1*theta3dot^2*sin(theta2 - theta3) - 18*l3*m2*theta3dot^2*sin(theta2 - theta3) - 4*l3*m3*theta3dot^2*sin(theta2 - theta3) + 6*l2*m2*theta2dot^2*sin(2*theta1 - 2*theta2) - 4*l2*m1*theta2dot^2*sin(2*theta2 - 2*theta3) + 4*l2*m3*theta2dot^2*sin(2*theta1 - 2*theta2) - 6*l2*m2*theta2dot^2*sin(2*theta2 - 2*theta3) - 4*l1*m1*theta1dot^2*sin(theta1 + theta2 - 2*theta3) - 6*l1*m2*theta1dot^2*sin(theta1 + theta2 - 2*theta3) - 6*l3*m2*theta3dot^2*sin(theta2 - 2*theta1 + theta3) - 4*l3*m3*theta3dot^2*sin(theta2 - 2*theta1 + theta3))/(2*l2*(2*m1 + 5*m2 + 2*m3 - 3*m2*cos(2*theta1 - 2*theta2) - 2*m1*cos(2*theta2 - 2*theta3) + m2*cos(2*theta1 - 2*theta3) - 2*m3*cos(2*theta1 - 2*theta2) - 3*m2*cos(2*theta2 - 2*theta3)))
      

--------------编辑-------- ----------

我无法以可用的形式获得解决方案，但我认为我可以利用 Oscar 方法的修改形式。 Oscar 的解决方案的问题在于，最终的解决方案太大而无法使用。

我也许可以简化 Oscar 方法中的 ps，然后在执行矩阵乘法之前将数值代入该表达式，这将给我一个数值解。此方法可用于odeint 调用的函数中，以对我的方程进行数值积分。挑战在于简化ps（长度超过180 万个字符）的表达式。现在，我将坚持使用 MATLAB 获得的解决方案。我没有测试我刚才建议的方法，但我发布了它以防其他人遇到这个问题。

Answer 1

我正在将我的评论扩展为答案。

在 sympy 1.7 中，有一个基于多项式域的矩阵的新实现（未记录且主要是实验性的）。虽然这并没有真正向用户宣传，但现在solve、linsolve 和solve_linear_system 使用它来求解线性方程组。这个想法是让矩阵的元素具有类似于 numpy 数组的 dtype 的元素，除了这里的 dtypes 被称为“域”并且基于诸如环和域之类的数学概念，例如ZZ[x, y] x 和 y 中的环多项式，具有整数系数。

由于浮点数以及 sin 和 cos 函数，您的方程组（还）不能用这个矩阵实现中的域表示，所以让我们用 Rationals 和 Symbols 替换它们：

In [2]: eqns = [nsimplify(eqn) for eqn in eqns]  # replace floats

In [3]: thetas = [theta1, theta2, theta3]
   ...: reps = {sin(t): s for s, t in zip(symbols('s1:4'), thetas)}
   ...: reps.update({cos(t): c for c, t in zip(symbols('c1:4'), thetas)})
   ...: eqns = [eq.subs(reps) for eq in eqns]    # replace sin/cos

我们现在可以为系统构造矩阵并将其转换为实验性的 DomainMatrix 类型：

In [5]: M, b = linear_eq_to_matrix(eqns, vars)

In [6]: from sympy.polys.domainmatrix import DomainMatrix

In [7]: dM = DomainMatrix.from_list_sympy(*M.shape, M.tolist())

In [8]: dM
Out[8]: DomainMatrix([[m1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0, 0, 0, 0], [0, m1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0], [0, 0, 0, m2, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1, 0, 0, 0], [0, 0, 0, 0, m2, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, m3, 0, 0, 0, 0, 0, 0, 1, 0, -1, 0], [0, 0, 0, 0, 0, 0, 0, m3, 0, 0, 0, 0, 0, 0, -1, 0, -1], [-1/2*c1*l1*m1, -1/2*s1*l1*m1, -1/12*l1**2*m1, 0, 0, 0, 0, 0, 0, 0, 0, c1*l1, -s1*l1, 0, 0, 0, 0], [0, 0, 0, -1/2*c2*l2*m2, -1/2*s2*l2*m2, -1/12*l2**2*m2, 0, 0, 0, 0, 0, 0, 0, c2*l2, -s2*l2, 0, 0], [0, 0, 0, 0, 0, 0, -1/2*c3*l3*m3, -1/2*s3*l3*m3, -1/12*l3**2*m3, 0, 0, 0, 0, 0, 0, c3*l3, s3*l3], [1, 0, -1/2*c1*l1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, -1/2*s1*l1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, -c1*l1, 1, 0, -1/2*c2*l2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, -s1*l1, 0, 1, -1/2*s2*l2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, -c1*l1, 0, 0, -c2*l2, 1, 0, -1/2*c3*l3, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, -s1*l1, 0, 0, -s2*l2, 0, 1, -1/2*s3*l3, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, c1*l1, 0, 0, c2*l2, 0, 0, c3*l3, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, s1*l1, 0, 0, s2*l2, 0, 0, s3*l3, 0, 0, 0, 0, 0, 0, 0, 0]], (17, 17), QQ[s1,s2,s3,c1,c2,c3,l1,l2,l3,m1,m2,m3])

In [9]: dM.domain
Out[9]: QQ[s1,s2,s3,c1,c2,c3,l1,l2,l3,m1,m2,m3]

我们看到，域是给定符号 l1, l2, ... 和符号 c1, s1 等中的有理数 QQ 上的多项式环，我们刚刚替换了 cos(theta1), sin(theta1), ...。

现在在求解的 1.7 实现中，将在相应的场（而不是环）上构建增广矩阵，并且系统将使用高斯消元法求解，仅使用基本的旋转（这里有改进的空间）。

由于这是一个方形系统，因此有一种可能更有效的替代方法。首先我们计算系统的特征多项式：

In [10]: %time p = dM.charpoly()
CPU times: user 1min 2s, sys: 1.56 s, total: 1min 4s
Wall time: 1min 8s

特征多项式的系数表达式很复杂，但基于 Berkowitz 算法相对较快： https://en.wikipedia.org/wiki/Samuelson%E2%80%93Berkowitz_algorithm

现在我们可以使用特征多项式和矩阵向量乘法来求解线性系统。为此，只需使用原始矩阵M 就很方便，因此我们不需要将域与b 统一。首先，我们将特征多项式的系数转换回正常的 sympy 表达式（来自环）：

In [15]: ps = [dM.domain.to_sympy(pi) for pi in p]

In [18]: sol = ps[0] * b

In [19]: for n in range(1, 17):
    ...:     sol = M*sol + ps[n]*b
    ...: 

In [20]: sol = sol / (-ps[17])  # Divide by the determinant

我可能在这方面犯了一个错误:)所以如果你打算使用这样的东西，也许值得在一个更简单的矩阵上检查这些计算。

尽管有一个重要的警告，但最终结果应该是一个解决方案。用符号 c 和 s 替换 sin(theta) 和 cos(theta) 忽略了这些在代数上依赖的事实，因为 c**2 + s**2 = 1。可以构建一个可以考虑这种代数依赖性的多项式域，并且它实际上可能更有效，因为会进行中间简化（例如，可以减少 c 的任何幂）。由于这种技术不涉及选择枢轴，唯一的风险是我们可能没有意识到行列式 ps[17] 何时为零。

上面显示的操作在我的计算机上总共需要大约 2 分钟或更短的时间。然而，最终结果以巨大的表达形式出现。如果我们在域中完成了sol 的最终构造，计算会更慢，但我们最终可能会得到更简单的最终结果。并非所有部分都在 1.7 中实现，但我可以轻松地向您展示该计算。