Does the code not has been optimized or my labtop too weak to run the code below?

Hi everybody,

I used the code below to solve Lagrange equation for Cantilever beam in dynamic case. I think the code has been written correctly, however when I ran this code, it stuck to line 80 (slover function) for few hours. I really appreciate if anyone can come up with any suggestion to modify the solver function or optimize the code.

syms M1 M2 M3 integer;
syms x1 x1d x1dd x2 x2d x2dd x3 x3d x3dd;
syms L L0 L1 L2 L3;
syms u1 u2 u3;
syms g m rho b h integer;
syms t integer;

Bn1 = [[cos(x1) -sin(x1) L1];[sin(x1) cos(x1) 0]; [0 0 1]];
Bn2 = [[cos(x2) -sin(x2) L2];[sin(x2) cos(x2) 0]; [0 0 1]];
Bn3 = [[cos(x3) -sin(x3) L3];[sin(x3) cos(x3) 0]; [0 0 1]];

B1 = Bn1;
B2 = Bn1.*Bn2;
B3 = Bn1.*Bn2.*Bn3;

Bd1 = diff(B1,x1)*x1d;
Bd2 = diff(B2,x1)*x1d + diff(B2,x2)*x2d;
Bd3 = diff(B3,x1)*x1d + diff(B3,x2)*x2d + diff(B3,x3)*x3d;

H1 = [[m*L1^2/12 0 0];[0 m*h^2/12 0]; [0 0 rho*b*h*L1]];
H2 = [[m*L2^2/12 0 0];[0 m*h^2/12 0]; [0 0 rho*b*h*L2]];
H3 = [[m*L3^2/12 0 0];[0 m*h^2/12 0]; [0 0 rho*b*h*L3]];

T1 = 1/2*trace(Bd1.*H1*Bd1.');
T1 = simplify(T1);

T2 = 1/2*trace(Bd2.*H2*Bd2.');
T2 = simplify(T2);

T3 = 1/2*trace(Bd3.*H3.*Bd3');
T3 = simplify(T3);
%Total Kinetic energy
T = T1 + T2 + T3;
%Total Potential energy
V = 0.5*5.12e-8*(x1^2 + x2^2 + x3^2);
%Energy dissipation due to rotational damping
D = 0.5*5.12e-17*(x1d^2 + x2d^2 + x3d^2);

Px1= u1;
Px2= u2;
Px3= u3;

pTpx1d = diff(T,x1d);
ddtpTpx1d = diff(pTpx1d,x1)*x1d+ ...
             diff(pTpx1d,x1d)*x1dd+ ...
             diff(pTpx1d,x2)*x2d + ...
             diff(pTpx1d,x2d)*x2dd+...
             diff(pTpx1d,x3)*x3d + ...
             diff(pTpx1d,x3d)*x3dd;
pTpx1 = diff(T,x1);
pVpx1 = diff(V,x1);
pDpx1 = diff(D,x1d);

pTpx2d = diff(T,x2d);
ddtpTpx2d = diff(pTpx2d,x1)*x1d+ ...
             diff(pTpx2d,x1d)*x1dd+ ...
             diff(pTpx2d,x2)*x2d + ...
             diff(pTpx2d,x2d)*x2dd+...
             diff(pTpx2d,x3)*x3d + ...
             diff(pTpx2d,x3d)*x3dd;
pTpx2 = diff(T,x2);
pVpx2 = diff(V,x2);
pDpx2 = diff(D,x2d);

pTpx3d = diff(T,x3d);
ddtpTpx3d = diff(pTpx3d,x1)*x1d+ ...
             diff(pTpx3d,x1d)*x1dd+ ...
             diff(pTpx3d,x2)*x2d + ...
             diff(pTpx3d,x2d)*x2dd+...
             diff(pTpx3d,x3)*x3d + ...
             diff(pTpx3d,x3d)*x3dd;
pTpx3 = diff(T,x3);
pVpx3 = diff(V,x3);
pDpx3 = diff(D,x3d);

eqx1 = simplify( ddtpTpx1d - pTpx1 + pVpx1 + pDpx1 - Px1);
eqx2 = simplify( ddtpTpx2d - pTpx2 + pVpx2 + pDpx2 - Px2);
eqx3 = simplify( ddtpTpx3d - pTpx3 + pVpx3 + pDpx3 - Px3);

Sol = solve(eqx1,eqx2,eqx3,'x1dd,x2dd,x3dd');
Sol.x1dd = simplify(Sol.x1dd);
Sol.x2dd = simplify(Sol.x2dd);
Sol.x3dd = simplify(Sol.x3dd);

syms y1 y2 y3 y4 y5 y6
fx1=subs(Sol.x1dd,{x1,x1d,x2,x2d,x3,x3d},{y1,y2,y3,y4,y5,y6})

Thanks, Tan

[Walter Roberson]

Which of the symbols can we assume are real-valued? Which can we assume are non-negative?

If you can add the assumption of real to all the symbols, then it simplifies the system enough to be able to solve. For example,

x1dd = (1/10141204801825835211973625643008)*(5070602400912917605986812821504*cos(x1)*((cos(x2)-1)*(cos(x2)+1)*(L3^2+h^2)*cos(x3)^4+((L2^2-L3^2)*cos(x2)^2-sin(x2)^2*(L3^2+h^2)*sin(x3)^2-L2^2+L3^2)*cos(x3)^2-(cos(x2)-1)*(cos(x2)+1)*(L2^2+h^2))*(L1^2+h^2)*m*x1d^2*sin(2*x1)+10141204801825835211973625643008*cos(x1)*(cos(x2)^2+sin(x2)^2-1)*m*sin(x1)*cos(x2)*(L3^2+h^2)^2*(cos(x1)*(x1d^2+x2d^2+x3d^2)*cos(x2)-2*x1d*x2d*sin(x1)*sin(x2))*cos(x3)^6-20282409603651670423947251286016*x3d*cos(x1)*cos(x2)*sin(x1)*m*sin(x3)*(cos(x2)^2+sin(x2)^2-1)*(L3^2+h^2)^2*(x1d*cos(x2)*sin(x1)+x2d*cos(x1)*sin(x2))*cos(x3)^5+10141204801825835211973625643008*(L3^2+h^2)*(cos(x1)^2*m*sin(x1)*((x1d^2+x2d^2+x3d^2)*(L3^2+h^2)*sin(x3)^2+(x1d^2+x2d^2)*h^2+(-x1d^2-x2d^2-x3d^2)*L3^2+2*L2^2*(x1d^2+x2d^2+(1/2)*x3d^2))*cos(x2)^4-2*((L3^2+h^2)*sin(x3)^2+h^2+2*L2^2-L3^2)*cos(x1)*x2d*m*sin(x1)^2*x1d*sin(x2)*cos(x2)^3+((sin(x2)+1)*m*sin(x1)*((x1d^2+x2d^2+x3d^2)*(L3^2+h^2)*sin(x3)^2+(x1d^2+x2d^2)*h^2+(-x1d^2-x2d^2-x3d^2)*L3^2+2*L2^2*(x1d^2+x2d^2+(1/2)*x3d^2))*(sin(x2)-1)*cos(x1)-12*u1+(23211375736600881/37778931862957161709568)*x1+(6230756230241793/10141204801825835211973625643008)*x1d)*cos(x1)*cos(x2)^2-2*sin(x1)*(((L3^2+h^2)*sin(x3)^2+h^2+2*L2^2-L3^2)*(sin(x2)+1)*x2d*m*sin(x1)*x1d*(sin(x2)-1)*cos(x1)-(23211375736600881/75557863725914323419136)*x2-(6230756230241793/20282409603651670423947251286016)*x2d+6*u2)*sin(x2)*cos(x2)+cos(x1)*(12*u1-(23211375736600881/37778931862957161709568)*x1-(6230756230241793/10141204801825835211973625643008)*x1d))*cos(x3)^4-20282409603651670423947251286016*sin(x3)*(cos(x2)^2+sin(x2)^2-1)*sin(x1)*(cos(x1)*((L3^2+h^2)*sin(x3)^2+L2^2-L3^2)*m*sin(x1)*x3d*x1d*cos(x2)^2+cos(x1)^2*((L3^2+h^2)*sin(x3)^2+L2^2-L3^2)*x2d*m*x3d*sin(x2)*cos(x2)-(23211375736600881/75557863725914323419136)*x3-(6230756230241793/20282409603651670423947251286016)*x3d+6*u3)*(L3^2+h^2)*cos(x3)^3+(10141204801825835211973625643008*cos(x1)^2*m*sin(x1)*(L2^2+h^2)*((x1d^2+x2d^2+x3d^2)*(L3^2+h^2)*sin(x3)^2+(-x1d^2-x2d^2-x3d^2)*h^2+(-2*x1d^2-2*x2d^2-x3d^2)*L3^2+L2^2*(x1d^2+x2d^2))*cos(x2)^4-20282409603651670423947251286016*cos(x1)*x2d*m*sin(x1)^2*((L3^2+h^2)*sin(x3)^2-h^2+L2^2-2*L3^2)*(L2^2+h^2)*x1d*sin(x2)*cos(x2)^3+10141204801825835211973625643008*cos(x1)*((sin(x2)+1)*m*sin(x1)*(L2^2+h^2)*((x1d^2+x2d^2+x3d^2)*(L3^2+h^2)*sin(x3)^2+(-x1d^2-x2d^2-x3d^2)*h^2+(-2*x1d^2-2*x2d^2-x3d^2)*L3^2+L2^2*(x1d^2+x2d^2))*(sin(x2)-1)*cos(x1)+(6230756230241793/10141204801825835211973625643008)*(L2-L3)*(L2+L3)*(-(40564819207303340847894502572032/2076918743413931)*u1+(2076918743413931127078912/2076918743413931)*x1+x1d))*cos(x2)^2-20282409603651670423947251286016*sin(x1)*((sin(x2)+1)*x2d*m*sin(x1)*((L3^2+h^2)*sin(x3)^2-h^2+L2^2-2*L3^2)*(L2^2+h^2)*x1d*(sin(x2)-1)*cos(x1)-(23211375736600881/75557863725914323419136)*(x2+(2076918743413931/2076918743413931127078912)*x2d-(151115727451828646838272/7737125245533627)*u2)*((L3^2+h^2)*sin(x3)^2+L2^2-L3^2))*sin(x2)*cos(x2)-6230756230241793*cos(x1)*(sin(x2)^2*(L3^2+h^2)*sin(x3)^2+L2^2-L3^2)*(-(40564819207303340847894502572032/2076918743413931)*u1+(2076918743413931127078912/2076918743413931)*x1+x1d))*cos(x3)^2-20282409603651670423947251286016*(x1d*x3d*cos(x1)*sin(x1)*m*(sin(x3)-1)*(sin(x3)+1)*(L3^2+h^2)*cos(x2)^2+x2d*x3d*cos(x1)^2*sin(x2)*m*(sin(x3)-1)*(sin(x3)+1)*(L3^2+h^2)*cos(x2)-(23211375736600881/75557863725914323419136)*x3-(6230756230241793/20282409603651670423947251286016)*x3d+6*u3)*sin(x3)*(cos(x2)^2+sin(x2)^2-1)*sin(x1)*(L2^2+h^2)*cos(x3)-10141204801825835211973625643008*(cos(x1)^2*sin(x1)*m*(x1d^2+x2d^2)*(L2^2+h^2)*cos(x2)^4-2*x1d*x2d*cos(x1)*sin(x1)^2*sin(x2)*m*(L2^2+h^2)*cos(x2)^3+cos(x1)*(sin(x1)*m*(sin(x2)-1)*(sin(x2)+1)*(x1d^2+x2d^2)*(L2^2+h^2)*cos(x1)-12*u1+(23211375736600881/37778931862957161709568)*x1+(6230756230241793/10141204801825835211973625643008)*x1d)*cos(x2)^2-2*sin(x1)*(x1d*x2d*sin(x1)*m*(sin(x2)-1)*(sin(x2)+1)*(L2^2+h^2)*cos(x1)-(23211375736600881/75557863725914323419136)*x2-(6230756230241793/20282409603651670423947251286016)*x2d+6*u2)*sin(x2)*cos(x2)+cos(x1)*(12*u1-(23211375736600881/37778931862957161709568)*x1-(6230756230241793/10141204801825835211973625643008)*x1d))*(L2^2+h^2))/(cos(x1)*m*(((cos(x1)^2-1)*cos(x2)^2-sin(x1)^2*sin(x2)^2-cos(x1)^2+1)*cos(x2)^2*(L3^2+h^2)^2*cos(x3)^6+(((2*L2^2-L3^2+h^2)*cos(x1)^2-(L3^2+h^2)*sin(x3)^2*sin(x1)^2-h^2-2*L2^2+L3^2)*cos(x2)^4+((-sin(x2)^2*(L3^2+h^2)*sin(x3)^2+L1^2-2*L2^2+L3^2)*cos(x1)^2+(-2*(sin(x2)^2-1/2)*(L3^2+h^2)*sin(x3)^2-sin(x2)^2*(2*L2^2-L3^2+h^2))*sin(x1)^2+sin(x2)^2*(L3^2+h^2)*sin(x3)^2-L1^2+2*L2^2-L3^2)*cos(x2)^2-(cos(x1)-1)*(cos(x1)+1)*(L1^2+h^2))*(L3^2+h^2)*cos(x3)^4+(-(L2^2+h^2)*((-L2^2+2*L3^2+h^2)*cos(x1)^2+(L3^2+h^2)*sin(x3)^2*sin(x1)^2-h^2+L2^2-2*L3^2)*cos(x2)^4+((-sin(x2)^2*(L3^2+h^2)*(L2^2+h^2)*sin(x3)^2+h^4+(L2^2+L3^2)*h^2+(-L1^2+2*L2^2)*L3^2+L1^2*L2^2-L2^4)*cos(x1)^2-2*(L2^2+h^2)*((sin(x2)^2-1/2)*(L3^2+h^2)*sin(x3)^2-(1/2)*sin(x2)^2*(-L2^2+2*L3^2+h^2))*sin(x1)^2+sin(x2)^2*(L3^2+h^2)*(L2^2+h^2)*sin(x3)^2-h^4+(-L2^2-L3^2)*h^2+(L1^2-2*L2^2)*L3^2-L1^2*L2^2+L2^4)*cos(x2)^2-(L1^2+h^2)*(sin(x2)^2*(L3^2+h^2)*sin(x3)^2+L2^2-L3^2)*(cos(x1)+1)*(cos(x1)-1))*cos(x3)^2-(L2^2+h^2)*((cos(x1)-1)*(cos(x1)+1)*(L2^2+h^2)*cos(x2)^4+((L1^2-L2^2)*cos(x1)^2-sin(x2)^2*(L2^2+h^2)*sin(x1)^2-L1^2+L2^2)*cos(x2)^2-(cos(x1)-1)*(cos(x1)+1)*(L1^2+h^2))))

... and if that is a meaningful answer to you then you are a better mathematician than I am ;-)