Non-linear differential equation system problem

Question

I have to prepare a system of differential equations to control, and so I need to solve the system for position, velocity, angle and angular velocity of both jaw and pitch angles. The equations are described by (11) (12) and (13) in the link: Two-wheeled balancing LEGO robot

So, my group mates used other softwares to compute the result, and they did get an aswer that is kinda close to what we expected, its just that the calculation fails before the whole pattern could be observed. So, I tried to get the results in Mathematica, but as it seems, the system is complex enough that one simple NDSolve won't help.

I have tried the code bellow:

{g, mb, mw, r, l, w, ip, ij, ixx, iw, jm, ra, kb, kt, bm} = {9.8, 1.5,
0.2, 0.1, 0.2, 1, 10, 10, 5, 5, 10, 10, 10, 10, 0.5}; 

u = HeavisideTheta[t];
x2[t] := D[x1[t], t];
x4[t] := D[x3[t], t];
x6[t] := D[x5[t], t];

system = {(mb + 2 mw + 2 (iw + jm)/r^2)*D[x2[t], t] + 
2/r*(bm + kt*kb/ra)*D[x1[t], t] + (mb*l - 2*jm/r)*D[x6[t], t] - 
2/r*(bm + kt*kb/ra)*D[x5[t], t] - mb*l*x5[t]*D[x5[t], t]^2 == 
kt*u/(r*ra), (2 (mw + (iw + jm)/r^2)*w^2 + ixx*x5[t]^2 + ij + 
mb*l*x5[t]^2)*D[x4[t], t] + 
2*((mb*l^2 + ixx - ij)*x5[t]*D[x5[t], t] + 
w^2/r^2*(bm + kt*kb/ra))*D[x3[t], t] == 0,
(mb*l^2 + ip + 2*jm)*D[x6[t], t] + (mb*l - 2*jm/r)*D[x2[t], t] + 
2*(bm + kt*kb/ra)*D[x3[t], t] - 
2/r*(bm + kt*kb/ra)*D[x1[t], t] - (mb*l^2 + ixx - ij)*x5[t]*
D[x3[t], t]^2 - mb*g*l*x5[t] == -kt*u/ra, x1[0] == 0, 
x2[0] == 0, x3[0] == 0, x4[0] == 0, x5[0] == 0, x6[0] == 0};

NDSolve[system, {x1, x2, x3, x4, x5, x6}, t]

Which then returns that this is not an ordinary differential equation. Our teacher has suggested that we use the transformation x'(t)=y(t) so that the equations were all first order differential equations. Any help is appreciated.

Answer 1

First of all the transformation suggested by your teacher should be part of your system. The derivatives are treated as separate variables. So your system will look like:

{g, mb, mw, r, l, w, ip, ij, ixx, iw, jm, ra, kb, kt, bm} = {9.8, 1.5,
    0.2, 0.1, 0.2, 1, 10, 10, 5, 5, 10, 10, 10, 10, 0.5};
u[t_] := HeavisideTheta[t];

system = {(mb + 2 mw + 2 (iw + jm)/r^2)*D[x2[t], t] + 
     2/r*(bm + kt*kb/ra)*D[x1[t], t] + (mb*l - 2*jm/r)*D[x6[t], t] - 
     2/r*(bm + kt*kb/ra)*D[x5[t], t] - mb*l*x5[t]*D[x5[t], t]^2 == 
    kt*u[t]/(r*ra), (2 (mw + (iw + jm)/r^2)*w^2 + ixx*x5[t]^2 + ij + 
        mb*l*x5[t]^2)*D[x4[t], t] + 
     2*((mb*l^2 + ixx - ij)*x5[t]*D[x5[t], t] + 
        w^2/r^2*(bm + kt*kb/ra))*D[x3[t], t] == 
    0, (mb*l^2 + ip + 2*jm)*D[x6[t], t] + (mb*l - 2*jm/r)*
      D[x2[t], t] + 2*(bm + kt*kb/ra)*D[x3[t], t] - 
     2/r*(bm + kt*kb/ra)*D[x1[t], t] - (mb*l^2 + ixx - ij)*x5[t]*
      D[x3[t], t]^2 - mb*g*l*x5[t] == -kt*u[t]/ra, 
   x2[t] == D[x1[t], t], x4[t] == D[x3[t], t], x6[t] == D[x5[t], t] , 
   x1[0] == 0, x2[0] == 0.0, x3[0] == 0, x4[0] == 0.0, x5[0] == 0, 
   x6[0] == 0.0};

NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 10}]

The above gives error because Mathematica tries to differentiate HevisideTheta function at $t=0$

Per @J.M. suggestion I replaced HevisideTheta with u[t_] := UnitStep[t]

And now the NDSolve returns solution, although it complains about singularity. Please check your equations.