# Solving a system of non linear differential equation of 6 DOF

Posted 10 years ago
28092 Views
|
3 Replies
|
2 Total Likes
|
 Hi forumI'm a beginner user in Mathematica, I'm trying to solve a non linear DE of 6 DOF using NDSolve, let me say something about this systemThe governing differential equation isM.ddq + V.dq + C = QQ is defined as vector of generalized forcesddq is defined as a vector of acceleration (ddx,ddy,ddz)dq is defined as a vector of velocity (dx,dy,dz)q is a vector of position (x, y , z)M and V are matrixsC a vector of component of weight as a fucntion of position of (x,y,z)NDsolve solve by x, y, z, rotx, roty, rotzQ has a a spring forces and damping forces of the suspension system, defined in a general way asFs=k*qFd=c*dqMy issue is that the damper forces changes as the velocity vector dq changes, soIf dq ( y ' ) > vcritthenFd=SIGN( y ' ) * ConstelseifFd = c*dqwhere vcrit is a constant value. The y ' term is got by NDSolve, but I don't know how to "modify" the Fd value once the NDsolve is running.Here is a piece of code: V0d1 = d1J.dq; dLd1 = -V0d1; yd1u = dd1/Norm[dd1]; fd1 = -cd1*(dLd1.yd1u); Fd1 = fd1*yd1u; F0 = Fs1 + Fs2 + Fd1 + Fd2 + Fd3 + Fd4 Q = {F0[[1]], F0[[2]], F0[[3]], M0[[1]], M0[[2]], M0[[3]],     M0[[4]]} /. {x -> x[t], y -> y[t],     z -> z[t], \[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][      t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t],     dx -> x'[t], dy -> y'[t], dz -> z'[t], d\[Alpha] -> \[Alpha]'[t],    d\[Beta] -> \[Beta]'[t], d\[Gamma] -> \[Gamma]'[t],     d\[Theta] -> \[Theta]'[t], ddx -> x''[t], ddy -> y''[t],     ddz -> z''[t], dd\[Alpha] -> \[Alpha]''[t],     dd\[Beta] -> \[Beta]''[t], dd\[Gamma] -> \[Gamma]''[t],     dd\[Theta] -> \[Theta]''[t]};ec1 = MD.ddq + MV.dq - MC /. {x -> x[t], y -> y[t],     z -> z[t], \[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][      t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t],     dx -> x'[t], dy -> y'[t], dz -> z'[t], d\[Alpha] -> \[Alpha]'[t],     d\[Beta] -> \[Beta]'[t], d\[Gamma] -> \[Gamma]'[t],     d\[Theta] -> \[Theta]'[t], ddx -> x''[t], ddy -> y''[t],     ddz -> z''[t], dd\[Alpha] -> \[Alpha]''[t],     dd\[Beta] -> \[Beta]''[t], dd\[Gamma] -> \[Gamma]''[t],     dd\[Theta] -> \[Theta]''[t]};MemoryConstrained[  sol1 = NDSolve[{ec1[[1]] == Q[[1]], ec1[[2]] == Q[[2]],      ec1[[3]] == Q[[3]], ec1[[4]] == Q[[4]], ec1[[5]] == Q[[5]],      ec1[[6]] == Q[[6]], ec1[[7]] == Q[[7]], x[0] == 0, y[0] == 0,      z[0] == 0, \[Alpha][0] == 0, \[Beta][0] == 0, \[Gamma][0] ==       0, \[Theta][0] == 0, x'[0] == 0, y'[0] == 0,      z'[0] == 0, \[Alpha]'[0] == 0, \[Beta]'[0] == 0, \[Gamma]'[0] ==       0, \[Theta]'[0] == 0}, {x, y,      z, \[Alpha], \[Beta], \[Gamma], \[Theta]}, {t, 0, tf},     MaxSteps -> \[Infinity]], Infinity];I'll be very thankful is someone can help me.Thanks
3 Replies
Sort By:
Posted 10 years ago
 I'm working through your large notebook, trying to see if I can understand enough to give you some suggestions.In your latest comment there appear to be some font or paste issues with lots of examples of \ -> \There are also appear to be a number of functions, like Rz4, Rz5, Rz6, Sz1, Sz2, Sz3, etc. that I assume you have defined outside the notebook somewhere.Something I see you repeatedly do in your notebook are things likeJacm = { {Id3[[1,1]], Id3[[1,2]], Id3[[1,3]], RrmR1[[1,1]], RrmR1[[1,2]], RrmR1[[1,3]], RrmR1[[1,4]]}, {Id3[[2,1]], Id3[[2,2]], Id3[[2,3]], RrmR1[[2,1]], RrmR1[[2,2]], RrmR1[[2,3]], RrmR1[[2,4]]}, {Id3[[3,1]], Id3[[3,2]], Id3[[3,3]], RrmR1[[3,1]], RrmR1[[3,2]], RrmR1[[3,3]], RrmR1[[3,4]]}};You might consider whether it would be more compact and perhaps a little less error prone to use this.{*only need define these once*)matrixJoin[m1_,m2_] := MapThread[Join,{m1,m2}];matrixJoin[m1_,m2_,m3_,m4_] := Join[matrixJoin[m1,m2],matrixJoin[m3,m4]];(*and then use matrixJoin as needed*)Jacm = matrixJoin[Id3, RrmR1];Mak1 = matrixJoin[m1*Id3,m1*Rr1R1,RrTR10,RJ1R1];Please test this carefully to make certain that it is doing the same thing as your lists of individual element operations.I'm trying to see ways to make use of more vector and matrix operations, rather than individual element operations,to try to make the code more compact and possibly easier to understand.But these don't address the real question that you have.
Posted 10 years ago
 ok, thanks Bill for your comments, here is the code   m1 = 31.8152;      m2 = 11.251952;   mm = 0.680389;   fr = 1000;   fr1 = 1000;   Fdamper = 80 ;(*N*)   vcrit = 130;   (*Data for body 1 - stator *)  ks2 = 9.500*fr1;(*N /mm*)  ks1 = ks2;  lo = 188/fr; (*m*)   xcg1 = -3.65426/fr;    ycg1 = -3.367/fr;  zcg1 = -0.724087/fr;  (*Data for body 2*)    xcg2 = 0.0061/fr;  ycg2 = -14.7063/fr;  zcg2 = -100.2917/fr;  (*Data for body m - imbalance mass *)    xcgm = 230.56/fr;  ycgm = -120.07/fr;  zcgm = 102.09/fr;  (*Data for spring 1 *)(*Spring location on body 1-right*)  x1s1 = 259.080/fr;  y1s1 = 219.4630/fr;  z1s1 = -26.7376/fr;  x0s1 = 306.4260/fr;  y0s1 = 449.5800/fr;(*Spring location on ground*)  z0s1 = -26.7376/fr;  (*Data for spring 2 - left*)  x1s2 = -259.080/fr;  y1s2 = 219.4630/fr;  z1s2 = -26.7376/fr;  x0s2 = -306.4260/fr;  y0s2 = 449.5800/fr;  z0s2 = -26.7376/fr;  (*Data for damper 1 on body 1 -FRD *)  x1d1 = 225.7326/fr;  y1d1 = -224.4966/fr ;(*Front Right Damper*)  z1d1 = 180/fr;  x0d1 = 300.9770/fr;  y0d1 = -431.2290/fr;  z0d1 = 210.6171/fr;  (*Data for damper 2 RRD*)  x1d2 = 225.7326/fr;  y1d2 = -240.776/fr;  z1d2 = -180.0/fr;  x0d2 = 300.9770/fr;  y0d2 = -447.5100/fr;  z0d2 = -149.3819/fr;  (*Data for damper 3 - FLD *)  x1d3 = -225.7326/fr;  y1d3 = -224.4966/fr;  z1d3 = 180.00/fr;  x0d3 = -300.977/fr;  y0d3 = -431.2290/fr;  z0d3 = 210.6181/fr;  (*Data for damper 4 - RLD*)  x1d4 = -225.7326/fr;  y1d4 = -240.7776/fr;  z1d4 = -180.00/fr;  x0d4 = -300.9770/fr;  y0d4 = -447.51/fr;  z0d4 = -149.3819/fr;  cd1 = -0.6154*fr1;(*N s /mm*)  (*Matrices de inercia centroidales*)  Ixx1 = 3506370/fr1^2;  Iyy1 = 1994860/fr1^2;  Izz1 = 2857680/fr1^2;  Ixy1 = 426122/fr1^2;  Ixz1 = 19905.8/fr1^2;  Iyz1 = 147236/fr1^2;  (*Matriz de inercia de rotor*)  Ixx2 = 506081.89/fr1^2;  Iyy2 = 503575.14/fr1^2;  Izz2 = 467128.95/fr1^2;  Ixy2 = -27.96/fr1^2;  Ixz2 = 8.70218/fr1^2;  Iyz2 = 4738.52/fr1^2;  (*Matriz de inercia de masa desbalanceo*)  Ixxm = 1773.24/fr1^2;  Iyym = 1505.01/fr1^2;  Izzm = 592.73/fr1^2;  Ixym = 214.25/fr1^2;  Ixzm = -30.11/fr1^2;  Iyzm = 130.80/fr1^2;  grav = {0, -9.810, 0};    ddq = {ddx, ddy, ddz, dd\[Alpha], dd\[Beta], dd\[Gamma], dd\[Theta]};    dq = {dx, dy, dz, d\[Alpha], d\[Beta], d\[Gamma], d\[Theta]};        r10 = {x, y, z};  Id3 = IdentityMatrix[3];  rcg1p = {xcg1, ycg1, zcg1};  R10 = Rz4[\[Alpha]].Rz5[\[Beta]].Rz6[\[Gamma]];  r0cg1 = r10 + R10.rcg1p;    (*Stator velocity measure on CS 0*)  Srcg1 = Sz1[xcg1] + Sz2[ycg1] + Sz3[zcg1];  R\[Phi]1 = {{Cos[\[Beta]]*Cos[\[Gamma]], Sin[\[Gamma]], 0,       0}, {-Cos[\[Beta]]*Sin[\[Gamma]], Cos[\[Gamma]], 0,       0}, {Sin[\[Beta]], 0, 1, 0}};  Rr1R1 = -R10.Srcg1.R\[Phi]1;  Jac1 = {{Id3[[1, 1]], Id3[[1, 2]], Id3[[1, 3]], Rr1R1[[1, 1]],       Rr1R1[[1, 2]], Rr1R1[[1, 3]], Rr1R1[[1, 4]]},     {Id3[[2, 1]], Id3[[2, 2]], Id3[[2, 3]], Rr1R1[[2, 1]],       Rr1R1[[2, 2]], Rr1R1[[2, 3]], Rr1R1[[2, 4]]},     {Id3[[3, 1]], Id3[[3, 2]], Id3[[3, 3]], Rr1R1[[3, 1]],       Rr1R1[[3, 2]], Rr1R1[[3, 3]], Rr1R1[[3, 4]]}};  V0cg1 = Jac1.dq;      (*Stator angular velocity measure on CS 1*)  ceros = ConstantArray[0, {3, 3}];   (*This matrix transforms the gc located at its cs to cs1 - data from \  cad*)  JB1 = {{ceros[[1, 1]], ceros[[1, 2]], ceros[[1, 3]],       R\[Phi]1[[1, 1]], R\[Phi]1[[1, 2]], R\[Phi]1[[1, 3]],       R\[Phi]1[[1, 4]]},     {ceros[[2, 1]], ceros[[2, 2]], ceros[[2, 3]], R\[Phi]1[[2, 1]],       R\[Phi]1[[2, 2]], R\[Phi]1[[2, 3]], R\[Phi]1[[2, 4]]},     {ceros[[3, 1]], ceros[[3, 2]], ceros[[3, 3]], R\[Phi]1[[3, 1]],       R\[Phi]1[[3, 2]], R\[Phi]1[[3, 3]], R\[Phi]1[[3, 4]]}};  \[Omega]01 = JB1.dq;    rcg2p = {xcg2, ycg2, zcg2};  R20 = Rz4[\[Alpha]].Rz5[\[Beta]].Rz6[\[Gamma]].Rz6[\[Theta]];  r0cg2 = r10 + R20.rcg2p;    (*Velocidad de centro de gravedad rotor*)  Srcg2 = Sz1[xcg2] + Sz2[ycg2] + Sz3[zcg2];  R\[Phi]2 = {{Cos[\[Beta]]*Cos[\[Gamma] + \[Theta]],       Sin[\[Gamma] + \[Theta]], 0,       0}, {-Cos[\[Beta]]*Sin[\[Gamma] + \[Theta]],       Cos[\[Gamma] + \[Theta]], 0, 0}, {Sin[\[Beta]], 0, 1, 1}};  Rr2R1 = -R20.Srcg2.R\[Phi]2;  Jac2 = {{Id3[[1, 1]], Id3[[1, 2]], Id3[[1, 3]], Rr2R1[[1, 1]],       Rr2R1[[1, 2]], Rr2R1[[1, 3]], Rr2R1[[1, 4]]},     {Id3[[2, 1]], Id3[[2, 2]], Id3[[2, 3]], Rr2R1[[2, 1]],       Rr2R1[[2, 2]], Rr2R1[[2, 3]], Rr2R1[[2, 4]]},     {Id3[[3, 1]], Id3[[3, 2]], Id3[[3, 3]], Rr2R1[[3, 1]],       Rr2R1[[3, 2]], Rr2R1[[3, 3]], Rr2R1[[3, 4]]}};  V0cg2 = Jac2.dq;    (*Velocidad angular del rotor medido en su cg*)  JB2 = {{ceros[[1, 1]], ceros[[1, 2]], ceros[[1, 3]],       R\[Phi]2[[1, 1]], R\[Phi]2[[1, 2]], R\[Phi]2[[1, 3]],       R\[Phi]2[[1, 4]]},     {ceros[[2, 1]], ceros[[2, 2]], ceros[[2, 3]], R\[Phi]2[[2, 1]],       R\[Phi]2[[2, 2]], R\[Phi]2[[2, 3]], R\[Phi]2[[2, 4]]},     {ceros[[3, 1]], ceros[[3, 2]], ceros[[3, 3]], R\[Phi]2[[3, 1]],       R\[Phi]2[[3, 2]], R\[Phi]2[[3, 3]], R\[Phi]2[[3, 4]]}};  \[Omega]02 = JB2.dq;    (***********************Cuerpo Masa Desbalanceo**************)  rcgmp = {xcgm, ycgm, zcgm};  r0cgm = r10 + R20.rcgmp;    (*Velocidad de centro de gravedad rotor*)  Srcgm = Sz1[xcgm] + Sz2[ycgm] + Sz3[zcgm];  RrmR1 = -R20.Srcgm.R\[Phi]2;  Jacm = {{Id3[[1, 1]], Id3[[1, 2]], Id3[[1, 3]], RrmR1[[1, 1]],       RrmR1[[1, 2]], RrmR1[[1, 3]], RrmR1[[1, 4]]},     {Id3[[2, 1]], Id3[[2, 2]], Id3[[2, 3]], RrmR1[[2, 1]],       RrmR1[[2, 2]], RrmR1[[2, 3]], RrmR1[[2, 4]]},     {Id3[[3, 1]], Id3[[3, 2]], Id3[[3, 3]], RrmR1[[3, 1]],       RrmR1[[3, 2]], RrmR1[[3, 3]], RrmR1[[3, 4]]}};  V0cgm = Jacm.dq;    (*Velocidad angular del rotor medido en su cg*)  \[Omega]2m = \[Omega]02;  (*RrTR01=-Transpose[R1].Transpose[Srcg1].Transpose[R01];*)  RrTR10 = m1*Transpose[Rr1R1];  J1 = Inercia[Ixx1, Iyy1, Izz1, Ixy1, Ixz1, Iyz1];  RJ1R1 = Transpose[R\[Phi]1].(m1*Transpose[Srcg1].Srcg1 + J1).R\[Phi]1;  Mak1 = {{m1*Id3[[1, 1]], m1*Id3[[1, 2]], m1*Id3[[1, 3]],       m1*Rr1R1[[1, 1]], m1*Rr1R1[[1, 2]], m1*Rr1R1[[1, 3]],       m1*Rr1R1[[1, 4]]},     {m1*Id3[[2, 1]], m1*Id3[[2, 2]], m1*Id3[[2, 3]], m1*Rr1R1[[2, 1]],       m1*Rr1R1[[2, 2]], m1*Rr1R1[[2, 3]], m1*Rr1R1[[2, 4]]},     {m1*Id3[[3, 1]], m1*Id3[[3, 2]], m1*Id3[[3, 3]], m1*Rr1R1[[3, 1]],       m1* Rr1R1[[3, 2]], m1*Rr1R1[[3, 3]], m1*Rr1R1[[3, 4]]},     {RrTR10[[1, 1]], RrTR10[[1, 2]], RrTR10[[1, 3]], RJ1R1[[1, 1]],       RJ1R1[[1, 2]], RJ1R1[[1, 3]], RJ1R1[[1, 4]]},     {RrTR10[[2, 1]], RrTR10[[2, 2]], RrTR10[[2, 3]], RJ1R1[[2, 1]],       RJ1R1[[2, 2]], RJ1R1[[2, 3]], RJ1R1[[2, 4]]},     {RrTR10[[3, 1]], RrTR10[[3, 2]], RrTR10[[3, 3]], RJ1R1[[3, 1]],       RJ1R1[[3, 2]], RJ1R1[[3, 3]], RJ1R1[[3, 4]]}, {RrTR10[[4, 1]],       RrTR10[[4, 2]], RrTR10[[4, 3]], RJ1R1[[4, 1]], RJ1R1[[4, 2]],       RJ1R1[[4, 3]], RJ1R1[[4, 4]]}};  (*K1=1/2*(dq.Mak1.dq);*)    (*Matriz de energiaa Cinetica cuerpo 2*)  RrTR20 = m2*Transpose[Rr2R1];  J2 = Inercia[Ixx2, Iyy2, Izz2, Ixy2, Ixz2, Iyz2];  RJ2R20 = Transpose[R\[Phi]2].(m2*Transpose[Srcg2].Srcg2 + J2).R\[Phi]2;  Mak2 = {{m2*Id3[[1, 1]], m2*Id3[[1, 2]], m2*Id3[[1, 3]],       m2*Rr2R1[[1, 1]], m2*Rr2R1[[1, 2]], m2*Rr2R1[[1, 3]],       m2*Rr2R1[[1, 4]]},     {m2*Id3[[2, 1]], m2*Id3[[2, 2]], m2*Id3[[2, 3]], m2*Rr2R1[[2, 1]],       m2*Rr2R1[[2, 2]], m2*Rr2R1[[2, 3]], m2*Rr2R1[[2, 4]]},     {m2*Id3[[3, 1]], m2*Id3[[3, 2]], m2*Id3[[3, 3]], m2*Rr2R1[[3, 1]],       m2* Rr2R1[[3, 2]], m2*Rr2R1[[3, 3]], m2*Rr2R1[[3, 4]]},     {RrTR20[[1, 1]], RrTR20[[1, 2]], RrTR20[[1, 3]], RJ2R20[[1, 1]],       RJ2R20[[1, 2]], RJ2R20[[1, 3]], RJ2R20[[1, 4]]},     {RrTR20[[2, 1]], RrTR20[[2, 2]], RrTR20[[2, 3]], RJ2R20[[2, 1]],       RJ2R20[[2, 2]], RJ2R20[[2, 3]], RJ2R20[[2, 4]]},     {RrTR20[[3, 1]], RrTR20[[3, 2]], RrTR20[[3, 3]], RJ2R20[[3, 1]],       RJ2R20[[3, 2]], RJ2R20[[3, 3]], RJ2R20[[3, 4]]}, {RrTR20[[4, 1]],       RrTR20[[4, 2]], RrTR20[[4, 3]], RJ2R20[[4, 1]], RJ2R20[[4, 2]],       RJ2R20[[4, 3]], RJ2R20[[4, 4]]}};  (*K2=1/2*(dq.Mak2.dq);*)    (*Matriz de energiaa Cinetica cuerpo masa de desbalanceo*)  RmrTR20 = mm*Transpose[RrmR1];  Jm = Inercia[Ixxm, Iyym, Izzm, Ixym, Ixzm, Iyzm];  RJmR20 = Transpose[R\[Phi]2].(mm*Transpose[Srcgm].Srcgm + Jm).R\[Phi]2;  Makm = {{mm*Id3[[1, 1]], mm*Id3[[1, 2]], mm*Id3[[1, 3]],       mm*RrmR1[[1, 1]], mm*RrmR1[[1, 2]], mm*RrmR1[[1, 3]],       mm*RrmR1[[1, 4]]},     {mm*Id3[[2, 1]], mm*Id3[[2, 2]], mm*Id3[[2, 3]], mm*RrmR1[[2, 1]],       mm*RrmR1[[2, 2]], mm*RrmR1[[2, 3]], mm*RrmR1[[2, 4]]},     {mm*Id3[[3, 1]], mm*Id3[[3, 2]], mm*Id3[[3, 3]], mm*RrmR1[[3, 1]],       mm* RrmR1[[3, 2]], mm*RrmR1[[3, 3]], mm*RrmR1[[3, 4]]},     {RmrTR20[[1, 1]], RmrTR20[[1, 2]], RmrTR20[[1, 3]],       RJmR20[[1, 1]], RJmR20[[1, 2]], RJmR20[[1, 3]], RJmR20[[1, 4]]},     {RmrTR20[[2, 1]], RmrTR20[[2, 2]], RmrTR20[[2, 3]],       RJmR20[[2, 1]], RJmR20[[2, 2]], RJmR20[[2, 3]], RJmR20[[2, 4]]},     {RmrTR20[[3, 1]], RmrTR20[[3, 2]], RmrTR20[[3, 3]], RJmR20[[3, 1]],       RJmR20[[3, 2]], RJmR20[[3, 3]],       RJmR20[[3, 4]]}, {RmrTR20[[4, 1]], RmrTR20[[4, 2]],       RmrTR20[[4, 3]], RJmR20[[4, 1]], RJmR20[[4, 2]], RJmR20[[4, 3]],       RJmR20[[4, 4]]}};  dM111 = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};    (*Termino dM/dt (M12)*)  dR10 = D[R10 /. {\[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][          t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t]},       t] /. {\[Alpha][t] -> \[Alpha], \[Beta][t] -> \[Beta], \[Gamma][        t] -> \[Gamma], \[Theta][t] -> \[Theta],   \!$$\*SuperscriptBox["\[Alpha]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Alpha],   \!$$\*SuperscriptBox["\[Beta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Beta],   \!$$\*SuperscriptBox["\[Gamma]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Gamma],   \!$$\*SuperscriptBox["\[Theta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Theta]};  dR\[Phi]1 =     D[R\[Phi]1 /. {\[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][          t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t]},       t] /. {\[Alpha][t] -> \[Alpha], \[Beta][t] -> \[Beta], \[Gamma][        t] -> \[Gamma], \[Theta][t] -> \[Theta],   \!$$\*SuperscriptBox["\[Alpha]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Alpha],   \!$$\*SuperscriptBox["\[Beta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Beta],   \!$$\*SuperscriptBox["\[Gamma]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Gamma],   \!$$\*SuperscriptBox["\[Theta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Theta]};  dM112 = -m1*(dR10.Srcg1.R\[Phi]1 + R10.Srcg1.dR\[Phi]1);    (*Termino dM/dt (M21)*)  dM121 = Transpose[dM112];    (*Termino dM/dt (M22)*)  dM122 = D[      RJ1R1 /. {\[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][          t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t]},       t] /. {\[Alpha][t] -> \[Alpha], \[Beta][t] -> \[Beta], \[Gamma][        t] -> \[Gamma], \[Theta][t] -> \[Theta],   \!$$\*SuperscriptBox["\[Alpha]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Alpha],   \!$$\*SuperscriptBox["\[Beta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Beta],   \!$$\*SuperscriptBox["\[Gamma]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Gamma],   \!$$\*SuperscriptBox["\[Theta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Theta]};    dMak1 = {{dM111[[1, 1]], dM111[[1, 2]], dM111[[1, 3]], dM112[[1, 1]],       dM112[[1, 2]], dM112[[1, 3]], dM112[[1, 4]]},      {dM111[[2, 1]], dM111[[2, 2]], dM111[[2, 3]], dM112[[2, 1]],       dM112[[2, 2]], dM112[[2, 3]], dM112[[2, 4]]},     {dM111[[3, 1]], dM111[[3, 2]], dM111[[3, 3]], dM112[[3, 1]],       dM112[[3, 2]], dM112[[3, 3]], dM112[[3, 4]]},     {dM121[[1, 1]], dM121[[1, 2]], dM121[[1, 3]], dM122[[1, 1]],       dM122[[1, 2]], dM122[[1, 3]], dM122[[1, 4]]},     {dM121[[2, 1]], dM121[[2, 2]], dM121[[2, 3]], dM122[[2, 1]],       dM122[[2, 2]], dM122[[2, 3]], dM122[[2, 4]]},     {dM121[[3, 1]], dM121[[3, 2]], dM121[[3, 3]], dM122[[3, 1]],       dM122[[3, 2]], dM122[[3, 3]], dM122[[3, 4]]}, {dM121[[4, 1]],       dM121[[4, 2]], dM121[[4, 3]], dM122[[4, 1]], dM122[[4, 2]],       dM122[[4, 3]], dM122[[4, 4]]}};    (*D1j=1/2*(Ddqcdqcj.(Mak1+Transpose[Mak1]));    V1j=1/2(Ddqcdqc1.(dMak1+Transpose[dMak1]));*)    (*Terminos de \[PartialD]L/\[PartialD]q*)  (*Termino \[PartialD]Mak1/\[PartialD]q (1,1)*)  dpMa1dpq1 = D[Mak1, x];  dpMa1dpq2 = D[Mak1, y];  dpMa1dpq3 = D[Mak1, z];  dpMa1dpq4 = D[Mak1, \[Alpha]];  dpMa1dpq5 = D[Mak1, \[Beta]];  dpMa1dpq6 = D[Mak1, \[Gamma]];  dpMa1dpq7 = D[Mak1, \[Theta]];  (*V1jp=1/2*(dq.dpMa1dpqj);*)    dpRcg1dpq1 = D[r0cg1, x];  dpRcg1dpq2 = D[r0cg1, y];  dpRcg1dpq3 = D[r0cg1, z];  dpRcg1dpq4 = D[r0cg1, \[Alpha]];  dpRcg1dpq5 = D[r0cg1, \[Beta]];  dpRcg1dpq6 = D[r0cg1, \[Gamma]];  dpRcg1dpq7 = D[r0cg1, \[Theta]];  dM211 = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};    (*Termino dM/dt (M12)*)  dR20 = D[R20 /. {\[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][          t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t]},       t] /. {\[Alpha][t] -> \[Alpha], \[Beta][t] -> \[Beta], \[Gamma][        t] -> \[Gamma], \[Theta][t] -> \[Theta],   \!$$\*SuperscriptBox["\[Alpha]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Alpha],   \!$$\*SuperscriptBox["\[Beta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Beta],   \!$$\*SuperscriptBox["\[Gamma]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Gamma],   \!$$\*SuperscriptBox["\[Theta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Theta]};  dR\[Phi]2 =     D[R\[Phi]2 /. {\[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][          t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t]},       t] /. {\[Alpha][t] -> \[Alpha], \[Beta][t] -> \[Beta], \[Gamma][        t] -> \[Gamma], \[Theta][t] -> \[Theta],   \!$$\*SuperscriptBox["\[Alpha]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Alpha],   \!$$\*SuperscriptBox["\[Beta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Beta],   \!$$\*SuperscriptBox["\[Gamma]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Gamma],   \!$$\*SuperscriptBox["\[Theta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Theta]};  dM212 = -m2*(dR20.Srcg2.R\[Phi]2 + R20.Srcg2.dR\[Phi]2);    (*Termino dM/dt (M21)*)  dM221 = Transpose[dM212];    (*Termino dM/dt (M22)*)  dM222 = D[      RJ2R20 /. {\[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][          t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t]},       t] /. {\[Alpha][t] -> \[Alpha], \[Beta][t] -> \[Beta], \[Gamma][        t] -> \[Gamma], \[Theta][t] -> \[Theta],   \!$$\*SuperscriptBox["\[Alpha]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Alpha],   \!$$\*SuperscriptBox["\[Beta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Beta],   \!$$\*SuperscriptBox["\[Gamma]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Gamma],   \!$$\*SuperscriptBox["\[Theta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Theta]};    dMak2 = {{dM211[[1, 1]], dM211[[1, 2]], dM211[[1, 3]], dM212[[1, 1]],       dM212[[1, 2]], dM212[[1, 3]], dM212[[1, 4]]},      {dM211[[2, 1]], dM211[[2, 2]], dM211[[2, 3]], dM212[[2, 1]],       dM212[[2, 2]], dM212[[2, 3]], dM212[[2, 4]]},     {dM211[[3, 1]], dM211[[3, 2]], dM211[[3, 3]], dM212[[3, 1]],       dM212[[3, 2]], dM212[[3, 3]], dM212[[3, 4]]},     {dM221[[1, 1]], dM221[[1, 2]], dM221[[1, 3]], dM222[[1, 1]],       dM222[[1, 2]], dM222[[1, 3]], dM222[[1, 4]]},     {dM221[[2, 1]], dM221[[2, 2]], dM221[[2, 3]], dM222[[2, 1]],       dM222[[2, 2]], dM222[[2, 3]], dM222[[2, 4]]},     {dM221[[3, 1]], dM221[[3, 2]], dM221[[3, 3]], dM222[[3, 1]],       dM222[[3, 2]], dM222[[3, 3]], dM222[[3, 4]]}, {dM221[[4, 1]],       dM221[[4, 2]], dM221[[4, 3]], dM222[[4, 1]], dM222[[4, 2]],       dM222[[4, 3]], dM222[[4, 4]]}};  (*Terminos de \[PartialD]L2/\[PartialD]q*)  (*Termino \[PartialD]Mak2/\[PartialD]q (1,1)*)  dpMa2dpq1 = D[Mak2, x];  dpMa2dpq2 = D[Mak2, y];  dpMa2dpq3 = D[Mak2, z];  dpMa2dpq4 = D[Mak2, \[Alpha]];  dpMa2dpq5 = D[Mak2, \[Beta]];  dpMa2dpq6 = D[Mak2, \[Gamma]];  dpMa2dpq7 = D[Mak2, \[Theta]];    dpRcg2dpq1 = D[r0cg2, x];  dpRcg2dpq2 = D[r0cg2, y];  dpRcg2dpq3 = D[r0cg2, z];  dpRcg2dpq4 = D[r0cg2, \[Alpha]];  dpRcg2dpq5 = D[r0cg2, \[Beta]];  dpRcg2dpq6 = D[r0cg2, \[Gamma]];  dpRcg2dpq7 = D[r0cg2, \[Theta]];    (*Rotor - cuerpo masa de desbalanceo*)  (*Termino dM/dt (M11)*)  dMm11 = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};    (*Termino dM/dt (M12)*)  dMm12 = -mm*(dR20.Srcgm.R\[Phi]2 + R20.Srcgm.dR\[Phi]2);    (*Termino dM/dt (M21)*)  dMm21 = Transpose[dMm12];    (*Termino dM/dt (M22)*)  dMm22 = D[      RJmR20 /. {\[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][          t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t]},       t] /. {\[Alpha][t] -> \[Alpha], \[Beta][t] -> \[Beta], \[Gamma][        t] -> \[Gamma], \[Theta][t] -> \[Theta],   \!$$\*SuperscriptBox["\[Alpha]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Alpha],   \!$$\*SuperscriptBox["\[Beta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Beta],   \!$$\*SuperscriptBox["\[Gamma]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Gamma],   \!$$\*SuperscriptBox["\[Theta]", "\[Prime]", MultilineFunction->None]$$[t] -> d\[Theta]};    dMakm = {{dMm11[[1, 1]], dMm11[[1, 2]], dMm11[[1, 3]], dMm12[[1, 1]],       dMm12[[1, 2]], dMm12[[1, 3]], dMm12[[1, 4]]},      {dMm11[[2, 1]], dMm11[[2, 2]], dMm11[[2, 3]], dMm12[[2, 1]],       dMm12[[2, 2]], dMm12[[2, 3]], dMm12[[2, 4]]},     {dMm11[[3, 1]], dMm11[[3, 2]], dMm11[[3, 3]], dMm12[[3, 1]],       dMm12[[3, 2]], dMm12[[3, 3]], dMm12[[3, 4]]},     {dMm21[[1, 1]], dMm21[[1, 2]], dMm21[[1, 3]], dMm22[[1, 1]],       dMm22[[1, 2]], dMm22[[1, 3]], dMm22[[1, 4]]},     {dMm21[[2, 1]], dMm21[[2, 2]], dMm21[[2, 3]], dMm22[[2, 1]],       dMm22[[2, 2]], dMm22[[2, 3]], dMm22[[2, 4]]},     {dMm21[[3, 1]], dMm21[[3, 2]], dMm21[[3, 3]], dMm22[[3, 1]],       dMm22[[3, 2]], dMm22[[3, 3]], dMm22[[3, 4]]}, {dMm21[[4, 1]],       dMm21[[4, 2]], dMm21[[4, 3]], dMm22[[4, 1]], dMm22[[4, 2]],       dMm22[[4, 3]], dMm22[[4, 4]]}};      (*Terminos de \[PartialD]L2/\[PartialD]q*)  (*Termino \[PartialD]Mak2/\[PartialD]q (1,1)*)    dpMamdpq1 = D[Makm, x];  dpMamdpq2 = D[Makm, y];  dpMamdpq3 = D[Makm, z];  dpMamdpq4 = D[Makm, \[Alpha]];  dpMamdpq5 = D[Makm, \[Beta]];  dpMamdpq6 = D[Makm, \[Gamma]];  dpMamdpq7 = D[Makm, \[Theta]];    dpRcgmdpq1 = D[r0cgm, x];  dpRcgmdpq2 = D[r0cgm, y];  dpRcgmdpq3 = D[r0cgm, z];  dpRcgmdpq4 = D[r0cgm, \[Alpha]];  dpRcgmdpq5 = D[r0cgm, \[Beta]];  dpRcgmdpq6 = D[r0cgm, \[Gamma]];  dpRcgmdpq7 = D[r0cgm, \[Theta]];    (*Fuerzas generaizadas*)    (*Fuerzas por resorte 1 -right front*)  (*Fuerzas por resorte 1 -right front*)  r0s1g = {x0s1, y0s1, z0s1}; (*medido inercialmente*)  (*vector pos, CS1,body 1,i-esim spring*)  r1s11p = {x1s1, y1s1, z1s1};   r0s11 = r10 + R10.r1s11p;  ds1 = r0s1g - r0s11;  ys1u = ds1/Norm[ds1];  Lo = lo; (*initial spring lenght*)  fs1 = ks1*(Norm[ds1] - Lo);  Fs1 = fs1*ys1u; (*Does the force needs to be defined in inertial CS?*)    (*Fuerzas por resortes*)  r0s2g = {x0s2, y0s2, z0s2}; (*medido inercialmente*)  (*vector pos, CS1,body 1,i-esim spring*)  r1s21p = {x1s2, y1s2, z1s2};   r0s21 = r10 + R10.r1s21p;  ds2 = r0s2g - r0s21;  ys2u = ds2/Norm[ds2];  Lo = lo; (*initial spring lenght*)  fs2 = ks2*(Norm[ds2] - Lo);  Fs2 = fs2*ys2u; (*Does the force needs to be defined in inertial CS?*)    (*Fuerzas por amortiguadores*)  (*Damper Right front damper*)  (*Fuerzas por amortiguadores*)  (*Damper Right front damper*)  rd11 = {x1d1, y1d1, z1d1};  rd1g = {x0d1, y0d1, z0d1};  dd1 = rd1g - (r10 + R10.rd11);    Sd11 = Sz1[x1d1] + Sz2[y1d1] + Sz3[z1d1];  Rd11R1 = -R10.Sd11.R\[Phi]1; (*write damper position as function of \  dqc*)  d1J = {{Id3[[1, 1]], Id3[[1, 2]], Id3[[1, 3]], Rd11R1[[1, 1]],      Rd11R1[[1, 2]], Rd11R1[[1, 3]], Rd11R1[[1, 4]]},    {Id3[[2, 1]], Id3[[2, 2]], Id3[[2, 3]], Rd11R1[[2, 1]],      Rd11R1[[2, 2]], Rd11R1[[2, 3]], Rd11R1[[2, 4]]},    {Id3[[3, 1]], Id3[[3, 2]], Id3[[3, 3]], Rd11R1[[3, 1]],      Rd11R1[[3, 2]], Rd11R1[[3, 3]], Rd11R1[[3, 4]]}};  V0d1 = d1J.dq;  dLd1 = -V0d1;  yd1u = dd1/Norm[dd1];  fd1 = -cd1*(dLd1.yd1u);  rd12 = {x1d2, y1d2, z1d2};  rd2g = {x0d2, y0d2, z0d2};  dd2 = rd2g - (r10 + R10.rd12);    Sd12 = Sz1[x1d2] + Sz2[y1d2] + Sz3[z1d2];  Rd12R1 = -R10.Sd12.R\[Phi]1; (*write damper position as function of \  dqc*)  d2J = {{Id3[[1, 1]], Id3[[1, 2]], Id3[[1, 3]], Rd12R1[[1, 1]],      Rd12R1[[1, 2]], Rd12R1[[1, 3]], Rd12R1[[1, 4]]},    {Id3[[2, 1]], Id3[[2, 2]], Id3[[2, 3]], Rd12R1[[2, 1]],      Rd12R1[[2, 2]], Rd12R1[[2, 3]], Rd12R1[[2, 4]]},    {Id3[[3, 1]], Id3[[3, 2]], Id3[[3, 3]], Rd12R1[[3, 1]],      Rd12R1[[3, 2]], Rd12R1[[3, 3]], Rd12R1[[3, 4]]}};  V0d2 = d2J.dq;  dLd2 = -V0d2;  yd2u = dd2/Norm[dd2];  fd2 = -cd1*(dLd2.yd2u);  (*Fd2=fd2*yd2u; fddddddddddddddddddddddddddddddddd*)    (*Damper Front Left damper*)  rd13 = {x1d3, y1d3, z1d3};  rd3g = {x0d3, y0d3, z0d3};  dd3 = rd3g - (r10 + R10.rd13);    Sd13 = Sz1[x1d3] + Sz2[y1d3] + Sz3[z1d3];  Rd13R1 = -R10.Sd13.R\[Phi]1; (*write damper position as function of \  dqc*)  d3J = {{Id3[[1, 1]], Id3[[1, 2]], Id3[[1, 3]], Rd13R1[[1, 1]],      Rd13R1[[1, 2]], Rd13R1[[1, 3]], Rd13R1[[1, 4]]},    {Id3[[2, 1]], Id3[[2, 2]], Id3[[2, 3]], Rd13R1[[2, 1]],      Rd13R1[[2, 2]], Rd13R1[[2, 3]], Rd13R1[[2, 4]]},    {Id3[[3, 1]], Id3[[3, 2]], Id3[[3, 3]], Rd13R1[[3, 1]],      Rd13R1[[3, 2]], Rd13R1[[3, 3]], Rd13R1[[3, 4]]}};  V0d3 = d3J.dq;  dLd3 = -V0d3;  yd3u = dd3/Norm[dd3];  fd3 = -cd1*(dLd3.yd3u);  (*Fd3=fd3*yd3u; fddddddddddddddddddddddddddddddddd*)(*Damper RearLeft damper*)rd14 = {x1d4, y1d4, z1d4};rd4g = {x0d4, y0d4, z0d4};dd4 = rd4g - (r10 + R10.rd14);Sd14 = Sz1[x1d4] + Sz2[y1d4] + Sz3[z1d4];Rd14R1 = -R10.Sd14.R\[Phi]1; (*write damper position as function of \dqc*)d4J = {{Id3[[1, 1]], Id3[[1, 2]], Id3[[1, 3]], Rd14R1[[1, 1]],    Rd14R1[[1, 2]], Rd14R1[[1, 3]], Rd14R1[[1, 4]]},  {Id3[[2, 1]], Id3[[2, 2]], Id3[[2, 3]], Rd14R1[[2, 1]],    Rd14R1[[2, 2]], Rd14R1[[2, 3]], Rd14R1[[2, 4]]},  {Id3[[3, 1]], Id3[[3, 2]], Id3[[3, 3]], Rd14R1[[3, 1]],    Rd14R1[[3, 2]], Rd14R1[[3, 3]], Rd14R1[[3, 4]]}};V0d4 = d4J.dq;dLd4 = -V0d4;yd4u = dd4/Norm[dd4];fd4 = -cd1*(dLd4.yd4u);(*Fd4=fd4*yd4u; fddddddddddddddddddddddddddddddddd*)MTR2\[Phi]R20 = Transpose[R\[Phi]2].Transpose[R20];M\[Omega] = {0, 0, -4.0}(*M0=Transpose[R10.R\[Phi]1].(MR10rs11.Fs1+MR10rs21.Fs2+MR10rd11.Fd1+\MR10rd12.Fd2+MR10rd13.Fd3+MR10rd14.Fd4)+MTR2\[Phi]R20.M\[Omega];*)(*Fexc=(R01.Transpose[R13].ku).M\[Omega];*)Ss1 = Sz1[x1s1] + Sz2[y1s1] + Sz3[z1s1];Ss2 = Sz1[x1s2] + Sz2[y1s2] + Sz3[z1s2];MR10rs11 = -Transpose[R\[Phi]1].Transpose[Ss1].Transpose[R10];MR10rs21 = -Transpose[R\[Phi]1].Transpose[Ss2].Transpose[R10];MR10rd11 = -Transpose[R\[Phi]1].Transpose[Sd11].Transpose[R10];MR10rd12 = -Transpose[R\[Phi]1].Transpose[Sd12].Transpose[R10];MR10rd13 = -Transpose[R\[Phi]1].Transpose[Sd13].Transpose[R10];MR10rd14 = -Transpose[R\[Phi]1].Transpose[Sd14].Transpose[R10];M0 = MR10rs11.Fs1 + MR10rs21.Fs2 + MR10rd11.Fd1 + MR10rd12.Fd2 +    MR10rd13.Fd3 + MR10rd14.Fd4 + MTR2\[Phi]R20.M\[Omega];(*Q={F0[[1]],F0[[2]],F0[[3]],M0[[1]],M0[[2]],M0[[3]],M0[[4]]}/.{x->x[\t],y->y[t],z->z[t],\[Alpha]->\[Alpha][t],\[Beta]->\[Beta][t],\[Gamma]-\>\[Gamma][t],\[Theta]->\[Theta][t],dx->x'[t],dy->y'[t],dz->z'[t],d\\[Alpha]->\[Alpha]'[t],d\[Beta]->\[Beta]'[t],d\[Gamma]->\[Gamma]'[t],\d\[Theta]->\[Theta]'[t],ddx->x''[t],ddy->y''[t],ddz->z''[t],dd\[Alpha]\->\[Alpha]''[t],dd\[Beta]->\[Beta]''[t],dd\[Gamma]->\[Gamma]''[t],dd\\[Theta]->\[Theta]''[t]};*)(*Ecuacion final*)Ddqcdqc1 = {1, 0, 0, 0, 0, 0, 0};Ddqcdqc2 = {0, 1, 0, 0, 0, 0, 0};Ddqcdqc3 = {0, 0, 1, 0, 0, 0, 0};Ddqcdqc4 = {0, 0, 0, 1, 0, 0, 0};Ddqcdqc5 = {0, 0, 0, 0, 1, 0, 0};Ddqcdqc6 = {0, 0, 0, 0, 0, 1, 0};Ddqcdqc7 = {0, 0, 0, 0, 0, 0, 1};D11 = Ddqcdqc1.(Mak1 + Transpose[Mak1]);D12 = Ddqcdqc2.(Mak1 + Transpose[Mak1]);D13 = Ddqcdqc3.(Mak1 + Transpose[Mak1]);D14 = Ddqcdqc4.(Mak1 + Transpose[Mak1]);D15 = Ddqcdqc5.(Mak1 + Transpose[Mak1]);D16 = Ddqcdqc6.(Mak1 + Transpose[Mak1]);D17 = Ddqcdqc7.(Mak1 + Transpose[Mak1]);D21 = Ddqcdqc1.(Mak2 + Transpose[Mak2]);D22 = Ddqcdqc2.(Mak2 + Transpose[Mak2]);D23 = Ddqcdqc3.(Mak2 + Transpose[Mak2]);D24 = Ddqcdqc4.(Mak2 + Transpose[Mak2]);D25 = Ddqcdqc5.(Mak2 + Transpose[Mak2]);D26 = Ddqcdqc6.(Mak2 + Transpose[Mak2]);D27 = Ddqcdqc7.(Mak2 + Transpose[Mak2]);Dm1 = Ddqcdqc1.(Makm + Transpose[Makm]);Dm2 = Ddqcdqc2.(Makm + Transpose[Makm]);Dm3 = Ddqcdqc3.(Makm + Transpose[Makm]);Dm4 = Ddqcdqc4.(Makm + Transpose[Makm]);Dm5 = Ddqcdqc5.(Makm + Transpose[Makm]);Dm6 = Ddqcdqc6.(Makm + Transpose[Makm]);Dm7 = Ddqcdqc7.(Makm + Transpose[Makm]);V11 = Ddqcdqc1.(dMak1 + Transpose[dMak1]);V12 = Ddqcdqc2.(dMak1 + Transpose[dMak1]);V13 = Ddqcdqc3.(dMak1 + Transpose[dMak1]);V14 = Ddqcdqc4.(dMak1 + Transpose[dMak1]);V15 = Ddqcdqc5.(dMak1 + Transpose[dMak1]);V16 = Ddqcdqc6.(dMak1 + Transpose[dMak1]);V17 = Ddqcdqc7.(dMak1 + Transpose[dMak1]);V21 = Ddqcdqc1.(dMak2 + Transpose[dMak2]);V22 = Ddqcdqc2.(dMak2 + Transpose[dMak2]);V23 = Ddqcdqc3.(dMak2 + Transpose[dMak2]);V24 = Ddqcdqc4.(dMak2 + Transpose[dMak2]);V25 = Ddqcdqc5.(dMak2 + Transpose[dMak2]);V26 = Ddqcdqc6.(dMak2 + Transpose[dMak2]);V27 = Ddqcdqc7.(dMak2 + Transpose[dMak2]);Vm1 = Ddqcdqc1.(dMakm + Transpose[dMakm]);Vm2 = Ddqcdqc2.(dMakm + Transpose[dMakm]);Vm3 = Ddqcdqc3.(dMakm + Transpose[dMakm]);Vm4 = Ddqcdqc4.(dMakm + Transpose[dMakm]);Vm5 = Ddqcdqc5.(dMakm + Transpose[dMakm]);Vm6 = Ddqcdqc6.(dMakm + Transpose[dMakm]);Vm7 = Ddqcdqc7.(dMakm + Transpose[dMakm]);V11p = dq.dpMa1dpq1;V12p = dq.dpMa1dpq2;V13p = dq.dpMa1dpq3;V14p = dq.dpMa1dpq4;V15p = dq.dpMa1dpq5;V16p = dq.dpMa1dpq6;V17p = dq.dpMa1dpq7;V21p = dq.dpMa2dpq1;V22p = dq.dpMa2dpq2;V23p = dq.dpMa2dpq3;V24p = dq.dpMa2dpq4;V25p = dq.dpMa2dpq5;V26p = dq.dpMa2dpq6;V27p = dq.dpMa2dpq7;Vm1p = dq.dpMa2dpq1;Vm2p = dq.dpMa2dpq2;Vm3p = dq.dpMa2dpq3;Vm4p = dq.dpMa2dpq4;Vm5p = dq.dpMa2dpq5;Vm6p = dq.dpMa2dpq6;Vm7p = dq.dpMa2dpq7;C11p = m1*(grav.dpRcg1dpq1);C12p = m1*(grav.dpRcg1dpq2);C13p = m1*(grav.dpRcg1dpq3);C14p = m1*(grav.dpRcg1dpq4);C15p = m1*(grav.dpRcg1dpq5);C16p = m1*(grav.dpRcg1dpq6);C17p = m1*(grav.dpRcg1dpq7);C21p = m2*(grav.dpRcg2dpq1);C22p = m2*(grav.dpRcg2dpq2);C23p = m2*(grav.dpRcg2dpq3);C24p = m2*(grav.dpRcg2dpq4);C25p = m2*(grav.dpRcg2dpq5);C26p = m2*(grav.dpRcg2dpq6);C27p = m2*(grav.dpRcg2dpq7);Cm1p = mm*(grav.dpRcgmdpq1);Cm2p = mm*(grav.dpRcgmdpq2);Cm3p = mm*(grav.dpRcgmdpq3);Cm4p = mm*(grav.dpRcgmdpq4);Cm5p = mm*(grav.dpRcgmdpq5);Cm6p = mm*(grav.dpRcgmdpq6);Cm7p = mm*(grav.dpRcgmdpq7);MD = 1/2*{D11 + D21 + Dm1, D12 + D22 + Dm2, D13 + D23 + Dm3,     D14 + D24 + Dm4, D15 + D25 + Dm5, D16 + D26 + Dm6,     D17 + D27 + Dm7};MV = 1/2*{V11 + V21 + Vm1 - (V11p + V21p + Vm1p),     V12 + V22 + Vm2 - (V12p + V22p + Vm2p),    V13 + V23 + Vm3 - (V13p + V23p + Vm3p),    V14 + V24 + Vm4 - (V14p + V24p + Vm4p),    V15 + V25 + Vm5 - (V15p + V25p + Vm5p),    V16 + V26 + Vm6 - (V16p + V26p + Vm6p),    V17 + V27 + Vm7 - (V17p + V27p + Vm7p)};MC = {C11p + C21p + Cm1p,   C12p + C22p + Cm2p,   C13p + C23p + Cm3p,   C14p + C24p + Cm4p,   C15p + C25p + Cm5p,   C16p + C26p + Cm6p,   C17p + C27p + Cm7p};ec1 = MD.ddq + MV.dq - MC /. {x -> x[t], y -> y[t],     z -> z[t], \[Alpha] -> \[Alpha][t], \[Beta] -> \[Beta][      t], \[Gamma] -> \[Gamma][t], \[Theta] -> \[Theta][t],     dx -> x'[t], dy -> y'[t], dz -> z'[t], d\[Alpha] -> \[Alpha]'[t],     d\[Beta] -> \[Beta]'[t], d\[Gamma] -> \[Gamma]'[t],     d\[Theta] -> \[Theta]'[t], ddx -> x''[t], ddy -> y''[t],     ddz -> z''[t], dd\[Alpha] -> \[Alpha]''[t],     dd\[Beta] -> \[Beta]''[t], dd\[Gamma] -> \[Gamma]''[t],     dd\[Theta] -> \[Theta]''[t]};tf = 15;MemoryConstrained[  sol1 = NDSolve[{ec1[[1]] == Q[[1]], ec1[[2]] == Q[[2]],      ec1[[3]] == Q[[3]], ec1[[4]] == Q[[4]], ec1[[5]] == Q[[5]],      ec1[[6]] == Q[[6]], ec1[[7]] == Q[[7]], x[0] == 0, y[0] == 0,      z[0] == 0, \[Alpha][0] == 0, \[Beta][0] == 0, \[Gamma][0] ==       0, \[Theta][0] == 0, x'[0] == 0, y'[0] == 0,      z'[0] == 0, \[Alpha]'[0] == 0, \[Beta]'[0] == 0, \[Gamma]'[0] ==       0, \[Theta]'[0] == 0}, {x, y,      z, \[Alpha], \[Beta], \[Gamma], \[Theta]}, {t, 0, tf},     MaxSteps -> \[Infinity]], Infinity];     On my understanding of how to programm a conditional I will do that,vcrit = 130;Fdam = 80;equaciones = {ec1[[1]] ==   If[(y' >= vcrit),    {Fd1 = Sign[y']*Fdam*yd1u;     Fd2 = Sign[y']*Fdam*yd2u;     Fd3 = Sign[y']*Fdam*yd3u;     Fd4 = Sign[y']*Fdam*yd4u;     F0 = Fs1 + Fs2 + Fd1 + Fd2 + Fd3 + Fd4;     Q = {F0[[1]], F0[[2]], F0[[3]], M0[[1]], M0[[2]], M0[[3]],         M0[[4]]} /. {x -> x, y -> y,         z -> z, \ -> \, \ -> \[          t], \ -> \, \ -> \,         dx -> x', dy -> y', dz -> z',         d\ -> \', d\ -> \',         d\ -> \', d\ -> \',         ddx -> x'', ddy -> y'', ddz -> z'',         dd\ -> \'', dd\ -> \'',         dd\ -> \'', dd\ -> \''};     },    {Fd1 = fd1*yd1u;     Fd2 = fd2*yd2u;     Fd3 = fd3*yd3u;     Fd4 = fd4*yd4u;     }    ]}But I know that in Mathematica there are some "methods" to stop the integration, so, as a begineer I don't how to do that.
Posted 10 years ago
 If you can share more of your problem then someone else might be able to experiment with it and find a solution. There are a handful of things I would try to see if I could make the code look more like the examples of differential equations that work, but, lacking any more information, I have no idea whether any of these would help you.