Why are there complex results for an inverted pendulum allowed to rock and slide with ODE113?

Hello everyone

I'm solving an inverted Pendulum Problem which is allowed to rock around the edges resp. allowed to slide. The Problem is defined in 3D.

While the solving of the Rocking Mode and the Sliding mode works well, there is a specific case where the combination of Rocking and Sliding combined is possible. I'm solving this using ODE113 but I always get complex results and I can't find out why. There are a few Squareroots in the ODE but they can't reach negative values, so there shouldn't be a reason there to get negative results.

There are 8 variables resp. 8 initial conditions which are all positive (but a small value around 1e-10).

Does anyone have an Idea where the problem is? Thank you in Advance!

The Code of the ODE is as followed:

function dy = cylinder_DGL_rockslide(t,y,time_r,dt,ug,vg,zg,Sys)
% CYLINDER_DGL evaluates the equation of motion for the rocking and rolling
% and sliding
% cylinder for given time instants t, excitation ug and vg and system
% properties R and H

%NEW excitation zg added

% Save geometry properties into new variables
r = Sys.r;
h = Sys.h;
g=9.81;
m = Sys.m;
mu = Sys.mu;


% Interpolate the excitation ug and vg for the time instant t
i = floor(t/dt);
if i < length(time_r)
    u = ug(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(ug(i+1)-ug(i));
    v = vg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(vg(i+1)-vg(i));

    %NEW ADDED z for vertical excitation
    w = zg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(zg(i+1)-zg(i));

else
    u = 0;
    v = 0;


    %NEW ADDED z for vertical excitation
    w = 0;
end


% Evaluate the equation of motion with help of the state space
% formulation
% State space formulation: row 1
dy(1,1)= y(2);

% State space formulation: row 2
dy(2,1)= 12.*m.^(-1).*((-10).*h.^2+(-9).*r.^2+6.*(h.^2+(-1).*r.^2).*cos(2.* ...
  y(1))+12.*h.*r.*sin(2.*y(1))).^(-1).*((1/24).*m.*((-12).*r.*(2.* ...
  h.*cos(y(1))+(-2).*h.*cos(2.*y(1))+3.*r.*sin(y(1)))+((-16).*h.^2+ ...
  15.*r.^2).*sin(2.*y(1))).*y(4).^2+m.*cos(y(3)).*( ...
  h.*cos(y(1))+r.*sin(y(1))).*u+m.*(h.*cos(y(1))+ ...
  r.*sin(y(1))).*sin(y(3)).*v+m.*(r.*cos(y(1))+( ...
  -1).*h.*sin(y(1))).*(g+w)+(-1).*m.*cos(y(3)).*( ...
  h.*cos(y(1))+r.*sin(y(1))).*(cos(y(3)).*(r.*cos(y(1))+(-1).*h.* ...
  sin(y(1))).*y(2).^2+(-2).*(h.*cos(y(1))+r.*sin(y(1) ...
  )).*sin(y(3)).*y(2).*y(4)+cos(y(3)) ...
  .*(r.*((-1)+cos(y(1)))+(-1).*h.*sin(y(1))).*y(4) ...
  .^2+u+mu.*y(6).*y(6).^2+y(8).^2).^(-1/2).*(g+w ...
  ))+(-1).*m.*(h.*cos(y(1))+r.*sin(y(1))).*sin(y(3)).*((r.*cos(y(1)) ...
  +(-1).*h.*sin(y(1))).*sin(y(3)).*y(2).^2+2.*cos(y(3).*(h.*cos(y(1))+r.*sin(y(1))).*y(2).*y(4)+(r.*((-1)+cos(y(1)))+(-1).*h.*sin(y(1))).*sin(y(3)).* ...
  y(4).^2+v+mu.*y(8) ...
  .*(y(6).^2+y(8).^2).^(-1/2).*(g+ ...
  w)));






% State space formulation: row 3
dy(3,1)= y(4);

% State space formulation: row 4
dy(4,1)= (-2).*(4.*h.^2+9.*r.^2+(4.*h.^2+(-3).*r.^2).*cos(y(1))).^(-1).*( ...
  y(6).^2+y(8).^2).^(-1/2).*((6.* ...
  r.^2+(4.*h.^2+(-3).*r.^2).*cos(y(1))).*cot((1/2).*y(1)).*( ...
  y(6).^2+y(8).^2).^(1/2).* ...
  y(2).*y(4)+6.*mu.*(r+h.*cot((1/2).*y(1))).*sin(y(3)).*y(6).*(g+w)+( ...
  -6).*mu.*cos(y(3)).*(r+h.*cot((1/2).*y(1))).*y(8).* ...
  (g+w));







% State space formulation: row 5
dy(5,1)= y(6);

% State space formulation: row 6
dy(6,1)= (1/16).*(4.*h.^2+9.*r.^2+(4.*h.^2+(-3).*r.^2).*cos(y(1))).^(-1).*( ...
  (-10).*h.^2+(-9).*r.^2+6.*(h.^2+(-1).*r.^2).*cos(2.*y(1))+12.*h.* ...
  r.*sin(2.*y(1))).^(-1).*(y(6).^2+y(8).^2).^(-1/2).*(4.*mu.*(448.*h.^4+1116.*h.^2.*r.^2+648.*r.^4+( ...
  352.*h.^4+(-90).*h.^2.*r.^2+(-450).*r.^4).*cos(y(1))+(-12).*(16.* ...
  h.^4+5.*h.^2.*r.^2+(-21).*r.^4).*cos(2.*y(1))+(-96).*h.^4.*cos(3.* ...
  y(1))+474.*h.^2.*r.^2.*cos(3.*y(1))+(-90).*r.^4.*cos(3.*y(1))+( ...
  -48).*h.^4.*cos(y(1)+(-2).*y(3))+69.*h.^2.*r.^2.*cos(y(1)+(-2).*y(3))+135.*r.^4.*cos(y(1)+(-2).*y(3))+48.*h.^4.*cos(3.*y(1)+(-2).*y(3))+(-237).*h.^2.*r.^2.*cos(3.*y(1)+(-2).*y(3))+45.*r.^4.*cos(3.* ...
  y(1)+(-2).*y(3))+96.*h.^4.*cos(2.*(y(1)+(-1).*y(3)))+30.*h.^2.* ...
  r.^2.*cos(2.*(y(1)+(-1).*y(3)))+(-126).*r.^4.*cos(2.*(y(1)+(-1).* ...
  y(3)))+(-192).*h.^4.*cos(2.*y(3))+(-300).*h.^2.*r.^2.*cos(2.*y(3)) ...
  +(-108).*r.^4.*cos(2.*y(3))+96.*h.^4.*cos(2.*(y(1)+y(3)))+30.* ...
  h.^2.*r.^2.*cos(2.*(y(1)+y(3)))+(-126).*r.^4.*cos(2.*(y(1)+y(3)))+ ...
  (-48).*h.^4.*cos(y(1)+2.*y(3))+69.*h.^2.*r.^2.*cos(y(1)+2.*y(3))+ ...
  135.*r.^4.*cos(y(1)+2.*y(3))+48.*h.^4.*cos(3.*y(1)+2.*y(3))+(-237) ...
  .*h.^2.*r.^2.*cos(3.*y(1)+2.*y(3))+45.*r.^4.*cos(3.*y(1)+2.*y(3))+ ...
  432.*h.^3.*r.*sin(y(1))+468.*h.*r.^3.*sin(y(1))+(-384).*h.^3.*r.* ...
  sin(2.*y(1))+(-504).*h.*r.^3.*sin(2.*y(1))+(-336).*h.^3.*r.*sin( ...
  3.*y(1))+324.*h.*r.^3.*sin(3.*y(1))+(-216).*h.^3.*r.*sin(y(1)+(-2) ...
  .*y(3))+(-234).*h.*r.^3.*sin(y(1)+(-2).*y(3))+168.*h.^3.*r.*sin( ...
  3.*y(1)+(-2).*y(3))+(-162).*h.*r.^3.*sin(3.*y(1)+(-2).*y(3))+192.* ...
  h.^3.*r.*sin(2.*(y(1)+(-1).*y(3)))+252.*h.*r.^3.*sin(2.*(y(1)+(-1) ...
  .*y(3)))+192.*h.^3.*r.*sin(2.*(y(1)+y(3)))+252.*h.*r.^3.*sin(2.*( ...
  y(1)+y(3)))+(-216).*h.^3.*r.*sin(y(1)+2.*y(3))+(-234).*h.*r.^3.* ...
  sin(y(1)+2.*y(3))+168.*h.^3.*r.*sin(3.*y(1)+2.*y(3))+(-162).*h.* ...
  r.^3.*sin(3.*y(1)+2.*y(3))).*y(6).*(g+w)+24.*mu.*((-32).*h.^4+(-50).*h.^2.*r.^2+(-18).*r.^4+((-16) ...
  .*h.^4+23.*h.^2.*r.^2+45.*r.^4).*cos(y(1))+2.*(16.*h.^4+5.*h.^2.* ...
  r.^2+(-21).*r.^4).*cos(2.*y(1))+16.*h.^4.*cos(3.*y(1))+(-79).* ...
  h.^2.*r.^2.*cos(3.*y(1))+15.*r.^4.*cos(3.*y(1))+(-72).*h.^3.*r.* ...
  sin(y(1))+(-78).*h.*r.^3.*sin(y(1))+64.*h.^3.*r.*sin(2.*y(1))+84.* ...
  h.*r.^3.*sin(2.*y(1))+56.*h.^3.*r.*sin(3.*y(1))+(-54).*h.*r.^3.* ...
  sin(3.*y(1))).*sin(2.*y(3)).*y(8).*(g+w)+2.*(y(6).^2+y(8).^2).^( ...
  1/2).*((-4).*(16.*h.^2+15.*r.^2).*cos(y(3)).*((-4).*h.^2.*r+3.* ...
  r.^3+(-2).*r.*(4.*h.^2+9.*r.^2).*cos(y(1))+((-4).*h.^2.*r+3.*r.^3) ...
  .*cos(2.*y(1))+8.*h.^3.*sin(y(1))+18.*h.*r.^2.*sin(y(1))+4.*h.^3.* ...
  sin(2.*y(1))+(-3).*h.*r.^2.*sin(2.*y(1))).*y(2).^2+ ...
  32.*r.*sin((1/2).*y(1)).*((-1).*(40.*h.^4+66.*h.^2.*r.^2+27.*r.^4) ...
  .*cos((1/2).*y(1))+3.*(4.*h.^4+(-13).*h.^2.*r.^2+(-3).*r.^4).*cos( ...
  (3/2).*y(1))+12.*h.^4.*cos((5/2).*y(1))+33.*h.^2.*r.^2.*cos((5/2) ...
  .*y(1))+(-9).*r.^4.*cos((5/2).*y(1))+60.*h.^3.*r.*sin((1/2).*y(1)) ...
  +54.*h.*r.^3.*sin((1/2).*y(1))+42.*h.^3.*r.*sin((3/2).*y(1))+6.* ...
  h.^3.*r.*sin((5/2).*y(1))+36.*h.*r.^3.*sin((5/2).*y(1))).*sin(y(3) ...
  ).*y(2).*y(4)+2.*(4.*h.^2+9.*r.^2+( ...
  4.*h.^2+(-3).*r.^2).*cos(y(1))).*(2.*cos(y(3)).*sin((1/2).*y(1)).* ...
  ((-2).*h.*(16.*h.^2+15.*r.^2).*cos((1/2).*y(1))+h.*(16.*h.^2+21.* ...
  r.^2).*cos((3/2).*y(1))+16.*h.^3.*cos((5/2).*y(1))+(-39).*h.* ...
  r.^2.*cos((5/2).*y(1))+(-40).*h.^2.*r.*sin((1/2).*y(1))+(-24).* ...
  r.^3.*sin((1/2).*y(1))+16.*h.^2.*r.*sin((3/2).*y(1))+21.*r.^3.* ...
  sin((3/2).*y(1))+40.*h.^2.*r.*sin((5/2).*y(1))+(-15).*r.^3.*sin(( ...
  5/2).*y(1))).*y(4).^2+(-4).*((-10).*h.^2+(-9).* ...
  r.^2+6.*(h.^2+(-1).*r.^2).*cos(2.*y(1))+12.*h.*r.*sin(2.*y(1))).* ...
  u+24.*cos(y(3)).*((-2).*h.*r.*cos(2.*y(1))+( ...
  h.^2+(-1).*r.^2).*sin(2.*y(1))).*(g+w))));



% State space formulation: row 7
dy(7,1)= y(8);

% State space formulation: row 8
dy(8,1)= (-1/16).*(4.*h.^2+9.*r.^2+(4.*h.^2+(-3).*r.^2).*cos(y(1))).^(-1).* ...
  ((-10).*h.^2+(-9).*r.^2+6.*(h.^2+(-1).*r.^2).*cos(2.*y(1))+12.*h.* ...
  r.*sin(2.*y(1))).^(-1).*(y(6).^2+y(8).^2).^(-1/2).*(24.*mu.*(32.*h.^4+50.*h.^2.*r.^2+18.*r.^4+(16.* ...
  h.^4+(-23).*h.^2.*r.^2+(-45).*r.^4).*cos(y(1))+(-2).*(16.*h.^4+5.* ...
  h.^2.*r.^2+(-21).*r.^4).*cos(2.*y(1))+(-16).*h.^4.*cos(3.*y(1))+ ...
  79.*h.^2.*r.^2.*cos(3.*y(1))+(-15).*r.^4.*cos(3.*y(1))+72.*h.^3.* ...
  r.*sin(y(1))+78.*h.*r.^3.*sin(y(1))+(-64).*h.^3.*r.*sin(2.*y(1))+( ...
  -84).*h.*r.^3.*sin(2.*y(1))+(-56).*h.^3.*r.*sin(3.*y(1))+54.*h.* ...
  r.^3.*sin(3.*y(1))).*sin(2.*y(3)).*y(6).*(g+ ...
  w)+4.*mu.*((-448).*h.^4+(-1116).*h.^2.*r.^2+( ...
  -648).*r.^4+((-352).*h.^4+90.*h.^2.*r.^2+450.*r.^4).*cos(y(1))+ ...
  12.*(16.*h.^4+5.*h.^2.*r.^2+(-21).*r.^4).*cos(2.*y(1))+96.*h.^4.* ...
  cos(3.*y(1))+(-474).*h.^2.*r.^2.*cos(3.*y(1))+90.*r.^4.*cos(3.*y(1))+(-48).*h.^4.*cos(y(1)+(-2).*y(3))+69.*h.^2.*r.^2.*cos(y(1)+( ...
  -2).*y(3))+135.*r.^4.*cos(y(1)+(-2).*y(3))+48.*h.^4.*cos(3.*y(1)+( ...
  -2).*y(3))+(-237).*h.^2.*r.^2.*cos(3.*y(1)+(-2).*y(3))+45.*r.^4.* ...
  cos(3.*y(1)+(-2).*y(3))+96.*h.^4.*cos(2.*(y(1)+(-1).*y(3)))+30.* ...
  h.^2.*r.^2.*cos(2.*(y(1)+(-1).*y(3)))+(-126).*r.^4.*cos(2.*(y(1)+( ...
  -1).*y(3)))+(-192).*h.^4.*cos(2.*y(3))+(-300).*h.^2.*r.^2.*cos(2.* ...
  y(3))+(-108).*r.^4.*cos(2.*y(3))+96.*h.^4.*cos(2.*(y(1)+y(3)))+ ...
  30.*h.^2.*r.^2.*cos(2.*(y(1)+y(3)))+(-126).*r.^4.*cos(2.*(y(1)+y(3)))+(-48).*h.^4.*cos(y(1)+2.*y(3))+69.*h.^2.*r.^2.*cos(y(1)+2.*y(3))+135.*r.^4.*cos(y(1)+2.*y(3))+48.*h.^4.*cos(3.*y(1)+2.*y(3))+( ...
  -237).*h.^2.*r.^2.*cos(3.*y(1)+2.*y(3))+45.*r.^4.*cos(3.*y(1)+2.* ...
  y(3))+(-432).*h.^3.*r.*sin(y(1))+(-468).*h.*r.^3.*sin(y(1))+384.* ...
  h.^3.*r.*sin(2.*y(1))+504.*h.*r.^3.*sin(2.*y(1))+336.*h.^3.*r.* ...
  sin(3.*y(1))+(-324).*h.*r.^3.*sin(3.*y(1))+(-216).*h.^3.*r.*sin(y(1)+(-2).*y(3))+(-234).*h.*r.^3.*sin(y(1)+(-2).*y(3))+168.*h.^3.* ...
  r.*sin(3.*y(1)+(-2).*y(3))+(-162).*h.*r.^3.*sin(3.*y(1)+(-2).*y(3) ...
  )+192.*h.^3.*r.*sin(2.*(y(1)+(-1).*y(3)))+252.*h.*r.^3.*sin(2.*(y(1)+(-1).*y(3)))+192.*h.^3.*r.*sin(2.*(y(1)+y(3)))+252.*h.*r.^3.* ...
  sin(2.*(y(1)+y(3)))+(-216).*h.^3.*r.*sin(y(1)+2.*y(3))+(-234).*h.* ...
  r.^3.*sin(y(1)+2.*y(3))+168.*h.^3.*r.*sin(3.*y(1)+2.*y(3))+(-162) ...
  .*h.*r.^3.*sin(3.*y(1)+2.*y(3))).*y(8).*(g+ ...
  w)+(y(6).^2+y(8) ...
  .^2).^(1/2).*(8.*(16.*h.^2+15.*r.^2).*((-4).*h.^2.*r+3.*r.^3+(-2) ...
  .*r.*(4.*h.^2+9.*r.^2).*cos(y(1))+((-4).*h.^2.*r+3.*r.^3).*cos(2.* ...
  y(1))+8.*h.^3.*sin(y(1))+18.*h.*r.^2.*sin(y(1))+4.*h.^3.*sin(2.*y(1))+(-3).*h.*r.^2.*sin(2.*y(1))).*sin(y(3)).*y(2) ...
  .^2+64.*r.*cos(y(3)).*sin((1/2).*y(1)).*((-1).*(40.*h.^4+66.* ...
  h.^2.*r.^2+27.*r.^4).*cos((1/2).*y(1))+3.*(4.*h.^4+(-13).*h.^2.* ...
  r.^2+(-3).*r.^4).*cos((3/2).*y(1))+12.*h.^4.*cos((5/2).*y(1))+33.* ...
  h.^2.*r.^2.*cos((5/2).*y(1))+(-9).*r.^4.*cos((5/2).*y(1))+60.* ...
  h.^3.*r.*sin((1/2).*y(1))+54.*h.*r.^3.*sin((1/2).*y(1))+42.*h.^3.* ...
  r.*sin((3/2).*y(1))+6.*h.^3.*r.*sin((5/2).*y(1))+36.*h.*r.^3.*sin( ...
  (5/2).*y(1))).*y(2).*y(4)+(-4).*(4.* ...
  h.^2+9.*r.^2+(4.*h.^2+(-3).*r.^2).*cos(y(1))).*(2.*sin((1/2).*y(1) ...
  ).*((-2).*h.*(16.*h.^2+15.*r.^2).*cos((1/2).*y(1))+h.*(16.*h.^2+ ...
  21.*r.^2).*cos((3/2).*y(1))+16.*h.^3.*cos((5/2).*y(1))+(-39).*h.* ...
  r.^2.*cos((5/2).*y(1))+(-40).*h.^2.*r.*sin((1/2).*y(1))+(-24).* ...
  r.^3.*sin((1/2).*y(1))+16.*h.^2.*r.*sin((3/2).*y(1))+21.*r.^3.* ...
  sin((3/2).*y(1))+40.*h.^2.*r.*sin((5/2).*y(1))+(-15).*r.^3.*sin(( ...
  5/2).*y(1))).*sin(y(3)).*y(4).^2+(-4).*((-10).* ...
  h.^2+(-9).*r.^2+6.*(h.^2+(-1).*r.^2).*cos(2.*y(1))+12.*h.*r.*sin( ...
  2.*y(1))).*v+24.*((-2).*h.*r.*cos(2.*y(1))+( ...
  h.^2+(-1).*r.^2).*sin(2.*y(1))).*sin(y(3)).*(g+w ...
  ))));
end 

[Torsten]

I suggest outputting "dy" in the course of the computation. Then you'll find out when and why you get complex results.

Answer 1

[Dan Dolan]

A number of the derivatives have factors to the 1/2 or -1/2 power. Whenver the terms inside become negative, imaginary results will appear. You either have to recheck the math for these derivatives or make sure to put absolute values on these roots.