Message Boards Message Boards

0
|
28471 Views
|
3 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Solving a system of non linear differential equation of 6 DOF

shair mendoza
Posted 11 years ago
Hi forum
I'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 system
The governing differential equation is
M.ddq + V.dq + C = Q

Q is defined as vector of generalized forces
ddq 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 matrixs
C a vector of component of weight as a fucntion of position of (x,y,z)

NDsolve solve by x, y, z, rotx, roty, rotz

Q has a a spring forces and damping forces of the suspension system, defined in a general way as
Fs=k*q
Fd=c*dq

My issue is that the damper forces changes as the velocity vector dq changes, so
If dq ( y ' ) > vcrit
then
Fd=SIGN( y ' ) * Const
elseif
Fd = c*dq

where 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
POSTED BY: shair mendoza
3 Replies
Sort By:
Bill Simpson
Bill Simpson
Posted 11 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 like
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]]}};
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 BY: Bill Simpson
shair mendoza
shair mendoza
Posted 11 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 BY: shair mendoza
Bill Simpson
Bill Simpson
Posted 11 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.
POSTED BY: Bill Simpson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract