Performance issues with solve in sympy

Question

First time posting. I am currently trying to solve 17 symbolic equations for 17 variables using sympy in Python. My equations are linear in the 17 variables. The issue that I am having is that the solve function is very slow - I have allowed my program to run for over two hours, and it still has not produced a result.

I have tried setting the flags "simplify" to false and "rational" to false, but my program still requires a very long time to run (it did not produce a result after several hours). I also tried using linsolve, and my program did not finish after I let it run overnight.

I also plugged my equations into MATLAB, converted them to matrix form, and solved them using MATLAB's backslash command. MATLAB gave me my result in a little over 2 minutes.

Does anyone have any advice on how to improve the performance of my Python code? I would like to be able to do all of my work in Python instead of having to rely on MATLAB's symbolic toolbox.

Python code:

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 code:

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 Output (This is only a portion of the output. The body of my post is limited to 30000 characters, and with the entire output posted, the post contains almost 67000 characters. I have done some checks, and this output is correct. It is also quite large - it is a 17 by 1 vector, with a newline separating each entry. I have pasted the first few entries)

(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)))
      

--------------------------Edit------------------------------

I have not been able to obtain a solution in a usable form, but I think that I might be able to make use of a modified form of Oscar's method. The issue with Oscar's solution is that the final solution is so massive that it is unusable.

I might be able to simplify the ps in Oscar's method and then substitute numerical values into this expression prior to performing matrix multiplication, which would give me a numerical solution. This method could be used in a function called by odeint to numerically integrate my equations. The challenge is in simplifying the expression for the ps (which is over 1.8 million characters long). For now, I will stick with the solution that I obtained using MATLAB. I have not tested the method that I just suggested, but I posted it in case anyone else runs into this issue.

Answer 1

I'm expanding my comment into an answer.

In sympy 1.7 there is a new (undocumented and largely experimental) implementation of matrices based on the polynomial domains. Although this isn't really advertised to users this is now used by solve, linsolve and solve_linear_system in order to solve systems of linear equations. The idea is to have matrices with elements that have something analogous to the dtype of a numpy array except that the dtypes here are called "domains" and are based on mathematical concepts like rings and fields e.g. ZZ[x, y] the ring polynomials in x and y with integer coefficients.

Your system of equations can not (yet) be represented by a domain in this matrix implementation because of the floats and the sin and cos functions so let's replace those with Rationals and 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

We can now construct the matrix for the system and convert it to the experimental DomainMatrix type:

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]

We see that the domain is a polynomial ring over the rationals QQ in the given symbols l1, l2, ... and also the symbols c1, s1 etc which we just substituted in for cos(theta1), sin(theta1), ....

Now in the 1.7 implementation of solve the augmented matrix would be constructed over the corresponding field (rather than ring) and the system would be solved using Gaussian elimination with only basic pivoting (there is room for improvement here).

Since this is a square system there is an alternative approach that is potentially more efficient. First we compute the characteristic polynomial of the system:

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

The expressions for the coefficients of the characteristic polynomial are complicated but this is relatively fast based on the Berkowitz algorithm: https://en.wikipedia.org/wiki/Samuelson%E2%80%93Berkowitz_algorithm

Now at this point we can solve the linear system using the characteristic polynomial and matrix-vector multiplications. For this it is convenient just to use the original matrix M so we don't need to unify the domain with b. First we convert the coefficients of the characteristic polynomial back to normal sympy expressions (from the ring):

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

I might have made a mistake in this :) so it's perhaps worth checking these calculations on a simpler matrix if you intend to use something like this.

The end result should be a solution although there is an important caveat. Replacing sin(theta) and cos(theta) with symbols c and s misses the fact that these are algebraically dependent because c**2 + s**2 = 1. It would be possible to build a polynomial domain that can take account of this algebraic dependence and it might actually be more efficient because intermediate simplification would be take place (e.g. any power of c could be reduced). Since this technique doesn't involve choosing pivots the only risk is that we might not realise when the determinant ps[17] is zero.

The operations shown above together take around 2 minutes or less on my computer. However the end result comes out in enormous expressions. Had we done the final construction of sol in the domains the calculation would have been slower but we might have ended up with a simpler final result. Not all the pieces are implemented in 1.7 for me to show that calculation easily to you though.