-7
$\begingroup$

I can't use my Mathematica 9 code in Mathematica 11 - it works in Mathematica 9, but when I execute it in Mathematica 11 it shows a lot of errors. Also, the code executes very slowly in Mathematica 9. If you can help me simplify my code and increase its efficiency, please help me.

Code:

g = 9.81;
a = 2;
b = 0;
R = (L1 + L2 + L3)/30;

L1 = 1;
L2 = 1;
L3 = 1;
L4 = 1;

LC1 = L1/2;
LC2 = L2/2;
LC3 = L3/2;
LC4 = L4/2;

(*Mass*)
m1 = 1;
m2 = 1;
m3 = 1;
m4 = 1;

(*Moment of inertia*)
I1 = 1;
I2 = 1;
I3 = 1;
I4 = 1;

(*Viscos damping factors*)
C1 = 1;
C2 = 1;
C3 = 1;
C4 = 1;


RDF = 1/2*C1*(θ1'[t])^2 + 1/2*C2*(θ2'[t])^2 + 
   1/2*C3*(θ3'[t])^2 + 1/2*C4*(θ4'[t])^2;

(*Input torques*)
T1 = 0;
T2 = 0;
T3 = 5*Sin[t];
T4 = 0;

xe[t_] := 
  a - L3*Cos[θ3[t]] - L4*Cos[θ4[t] - θ3[t]];
ye[t_] := -L3*Sin[θ3[t]] + L4*Sin[θ4[t] - θ3[t]];

θ2[t_] := 
  ArcCos[((xe[t])^2 + (ye[t] - b)^2 - L1^2 - L2^2)/(2 L1*L2)];

θ1[t_] := 
  ArcTan[(L1 + L2*Cos[θ2[t]])*(-xe[t]) + 
    L2*Sin[θ2[t]]*(-ye[t] - b), -L2*
     Sin[θ2[t]]*(-xe[t]) + (L1 + 
       L2*Cos[θ2[t]])*(-ye[t] - b)];

(*Center of gravity positions*)
XG1 = -LC1*Cos[θ1[t]];
YG1 = b - LC1*Sin[θ1[t]];
XG2 = -L1*Cos[θ1[t]] - LC2*Cos[θ1[t] + θ2[t]];
YG2 = b - L1*Sin[θ1[t]] - LC2*Sin[θ1[t] + θ2[t]];
XG3 = a - L3*Cos[θ3[t]];
YG3 = -LC3*Sin[θ3[t]];
XG4 = a - L3*Cos[θ3[t]] - LC4*Cos[θ4[t] - θ3[t]];
YG4 = -LC3*Sin[θ3[t]] + LC4*Sin[θ4[t] - θ3[t]];

ω1 = D[θ1[t], t];
ω2 = D[θ1[t] + θ2[t], t];
ω3 = D[θ3[t], t];
ω3 = D[θ4[t] - θ3[t], t];

T = 1/2*m1*((D[XG1, t])^2 + (D[YG1, t])^2) + 1/2*I1*ω1^2 + 
   1/2*m2*((D[XG2, t])^2 + (D[YG2, t])^2) + 1/2*I2*ω2^2 + 
   1/2*m3*((D[XG3, t])^2 + (D[YG3, t])^2) + 1/2*I3*ω3^2 + 
   1/2*m4*((D[XG4, t])^2 + (D[YG4, t])^2) + 1/2*I4*ω4^2;
V = m1*g*YG1 + m2*g*YG2 + m3*g*YG3 + m4*g*YG4;

(*Lagrangian*)
L = T - V;

(*Equations*)
(*Energy equations*)
eq1 = D[D[L, θ3'[t]], t] - D[L, θ3[t]] + 
    D[RDF, θ3'[t]] == T3;
eq2 = D[D[L, θ4'[t]], t] - D[L, θ4[t]] + 
    D[RDF, θ4'[t]] == T4;


(*eq1=θ1'[t]+6θ2'[t]+θ1''[t]\[Equal]0;
eq2=θ1'[t]-θ2'[t]+θ2''[t]\[Equal]0;*)


ics = {θ3[0] == Pi/2, θ3'[0] == 0, θ4[0] == 
    Pi/2, θ4'[0] == 0};
eqs = Join[{eq1, eq2}, ics];

s = NDSolve[eqs, {θ3, θ4}, {t, 0, 50}, 
   MaxSteps -> 1000000000];
$\endgroup$
4
  • 3
    $\begingroup$ Have you tried to debug it step by step and figure out when and why specific errors are printed? "fix my code and improve it later" - this site is not a free debugging service $\endgroup$
    – Kuba
    Feb 6, 2017 at 11:13
  • 1
    $\begingroup$ Did you really intend:ω3 = D[θ3[t], t]; ω3 = D[θ4[t] - θ3[t], t];? or should it be ω3 = D[θ3[t], t]; ω4 = D[θ4[t] - θ3[t], t]; otherwise ω4 is undefined $\endgroup$
    – Young
    Feb 6, 2017 at 22:07
  • $\begingroup$ @Farshid are you still waiting for an answer? $\endgroup$
    – zhk
    Feb 26, 2017 at 4:23
  • $\begingroup$ your answer was correct and I executed code easily. thanks for your help. $\endgroup$
    – Farshid
    Feb 28, 2017 at 19:32

1 Answer 1

8
$\begingroup$

The answer to your question is in the error,

NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

s = NDSolve[eqs, {θ3, θ4}, {t, 0, 50},MaxSteps -> 1000000000,
   Method -> {"EquationSimplification" -> "Residual"}];
Plot[{θ3[t], θ4[t]} /. s, {t, 0, 50}]

enter image description here

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.