在 sympy 中解决性能问题
Performance issues with solve in sympy
第一次posting。我目前正在尝试使用 Python 中的 sympy 求解 17 个变量的 17 个符号方程。我的方程式在 17 个变量中是线性的。我遇到的问题是解决函数非常慢 - 我已经让我的程序 运行 两个多小时了，但它仍然没有产生结果。
我已经尝试将标志“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 输出（这只是输出的一部分。我的 post 的正文限制为 30000 个字符，而整个输出 posted，post包含近 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 获得的解决方案。我还没有测试我刚刚建议的方法，但是我 post 编辑了它以防其他人 运行 遇到这个问题。
我正在将我的评论扩展为答案。
在 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, ... 中的有理数 QQ 以及我们刚刚替换为 [=24 的符号 c1, s1 等的多项式环=].
现在在 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 统一。首先我们将特征多项式的系数转换回正常的 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 中实现，但我可以轻松地向您展示该计算。
第一次posting。我目前正在尝试使用 Python 中的 sympy 求解 17 个变量的 17 个符号方程。我的方程式在 17 个变量中是线性的。我遇到的问题是解决函数非常慢 - 我已经让我的程序 运行 两个多小时了，但它仍然没有产生结果。
我已经尝试将标志“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 输出（这只是输出的一部分。我的 post 的正文限制为 30000 个字符，而整个输出 posted，post包含近 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 获得的解决方案。我还没有测试我刚刚建议的方法，但是我 post 编辑了它以防其他人 运行 遇到这个问题。
我正在将我的评论扩展为答案。
在 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, ... 中的有理数 QQ 以及我们刚刚替换为 [=24 的符号 c1, s1 等的多项式环=].
现在在 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 统一。首先我们将特征多项式的系数转换回正常的 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 中实现，但我可以轻松地向您展示该计算。