使用Symphy解决的性能问题

Question

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

我已经尝试将标志“simplify”设置为false，将“rational”设置为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乘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，然后在执行矩阵乘法之前将数值替换到这个表达式中，这将给我一个数值解对我的方程进行数值积分。挑战在于简化ps（长度超过180万个字符）的表达式。目前，我将坚持使用MATLAB获得的解决方案。我没有测试我刚才建议的方法，但我发布了它，以防其他人遇到此问题

Answer 1

我将我的评论扩展为一个答案

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

在这个矩阵实现中，由于浮点数和sin和cos函数，你的方程组还不能用域来表示，所以让我们用有理数和符号来代替它们：

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

我们现在可以为系统构建矩阵，并将其转换为实验域矩阵类型：

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), ...

现在在solve的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统一起来。首先，我们将特征多项式的系数转换回正常的辛表达式（从环）：

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

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

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

上面显示的操作在我的计算机上总共需要大约2分钟或更少的时间。然而，最终的结果以巨大的表达出来。如果我们在域中完成了sol的最终构造，那么计算会更慢，但我们可能会得到更简单的最终结果。但并不是所有的部分都在1.7中实现，我可以很容易地向您展示这个计算