第一次发布。我目前正在尝试在Python中使用sympy解决17个变量的17个符号方程式。我的方程在17个变量中是线性的。我遇到的问题是Solve函数的运行速度很慢-我已允许程序运行两个小时以上，但仍未产生结果。
我尝试将标志“简化”设置为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解决方案的问题在于最终解决方案规模巨大，无法使用。
我也许可以简化psOscar的方法，然后在执行矩阵乘法之前将数值替换为该表达式，这将为我提供数值解决方案。此方法可用于调用的函数中，odeint以对我的方程式进行数值积分。挑战在于简化表达式ps（长度超过180万个字符）。现在，我将坚持使用MATLAB获得的解决方案。我没有测试我刚刚建议的方法，但是我发布了它，以防其他人遇到此问题。
我将评论扩展为答案。
在sympy 1.7中，有一个基于多项式域的矩阵的新实现（未记录，且很大程度上是实验性的）。尽管并没有真正向用户宣传此内容，但现在已由solve，linsolve并且solve_linear_system为了解决线性方程组而使用了它。这个想法是让矩阵中的元素具有类似于dtypenumpy数组的元素的元素，只是这里的dtypes被称为“域”，并且基于诸如环和域之类的数学概念，例如整数系数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
现在，我们可以为系统构造矩阵并将其转换为实验性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]
我们看到，域是QQ给定符号l1, l2, ...以及c1, s1我们刚刚替换的符号等上有理数的多项式环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。首先，我们将特征多项式的系数转换回正常的对称表达式（从环）：
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中实现，但我还是可以轻松地向您展示该计算。
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句