ODE45 producing NaN values for a set of differential equations

Hi Guys

Please refer to the below code.

When I try and excute the function using ODE45 solver all the values produced are NaN.

What is the source of the issue? I have checked the code numerous times and there are no mistakes in equations.

clc;
clear;
r0 = 1;
L = 2.5;

k = 353160*2;
m1 = 3500/2;
m2 = 3500/2;


 SolverOptions = odeset('RelTol',1e-5,'AbsTol',1e-5);
 [t,y] = ode45(@Sspace,[0:0.001:10],[0 0 0 0 0 0],SolverOptions,r0,L,k,m1,m2);
 y

function xp = Sspace(t,x,r0,L,k,m1,m2)

% r2 = x(5);
% r2_dot = x(6);
% theta1 = x(1);
% theta1_dot = x(2);
% theta2 = x(3);
% theta2_dot = x(4);

one = sqrt(L+x(5)*cos(x(1))-x(5)*cos(x(3)));
two = sqrt(L-x(5)*cos(x(1))+x(5)*cos(x(3)));

r1 = (L+one*two+x(5)*sin(x(3)))/sin(x(1));

one_dot = x(6)*cos(x(1))-x(5)*x(2)*sin(x(1))-x(6)*cos(x(3))+x(5)*x(4)*sin(x(3));
two_dot = -x(6)*cos(x(1))+x(6)*x(2)*sin(x(1))+x(6)*cos(x(3))-x(5)*x(4)*sin(x(3));

R1_dot = one_dot*two/(2*one)+two_dot*one/(2*two);

r1_dot = ((R1_dot+x(6)*sin(x(3))+x(5)*x(4)*cos(x(3)))*sin(x(1))-(x(2)*cos(x(1))*(L+one*two+2*x(5)*sin(x(3)))))/sin(x(1))^2;

theta1_dotdot = (-2*m1*r1*r1_dot*x(2)-m1*9.81*r1*sin(x(1)))/(m1*r1^2);

theta2_dotdot = (-2*m2*x(5)*x(6)*x(4)-m2*9.81*x(5)*sin(x(3)))/(m2*x(5)^2);

r2_dotdot = (m2*x(5)*x(4)^2-k*(x(5)-r0)+m2*9.81*cos(x(3)))/m2;

xp = [x(2);
      theta2_dotdot;
      x(4);
      theta2_dotdot;
      x(6);
      r2_dotdot];
  

end

Answer 1

[Cris LaPierre]

Your numerator and denominators are both 0. And 0/0 = NaN

x=[0 0 0 0 0 0];
r0 = 1;
L = 2.5;
k = 353160*2;
m1 = 3500/2;
m2 = 3500/2;

one = sqrt(L+x(5)*cos(x(1))-x(5)*cos(x(3)));
two = sqrt(L-x(5)*cos(x(1))+x(5)*cos(x(3)));
r1 = (L+one*two+x(5)*sin(x(3)))/sin(x(1));
one_dot = x(6)*cos(x(1))-x(5)*x(2)*sin(x(1))-x(6)*cos(x(3))+x(5)*x(4)*sin(x(3));
two_dot = -x(6)*cos(x(1))+x(6)*x(2)*sin(x(1))+x(6)*cos(x(3))-x(5)*x(4)*sin(x(3));

R1_dot = one_dot*two/(2*one)+two_dot*one/(2*two);
r1_dot = ((R1_dot+x(6)*sin(x(3))+x(5)*x(4)*cos(x(3)))*sin(x(1))-(x(2)*cos(x(1))*(L+one*two+2*x(5)*sin(x(3)))))/sin(x(1))^2

r1_dot = NaN

% numerator
((R1_dot+x(6)*sin(x(3))+x(5)*x(4)*cos(x(3)))*sin(x(1))-(x(2)*cos(x(1))*(L+one*two+2*x(5)*sin(x(3)))))

ans = 0

% denominator
((R1_dot+x(6)*sin(x(3))+x(5)*x(4)*cos(x(3)))*sin(x(1))-(x(2)*cos(x(1))*(L+one*two+2*x(5)*sin(x(3)))))

ans = 0

% result of 0/0
0/0

ans = NaN

[Mishal Mohanlal]

Thank you kindly

Answer 2

[Paul]

One problem appears to be that the initial condition of x(5) is 0, which results in a divide by zero in the computation of theta2_dotdot.