ODE45 producing NaN values for a set of differential equations

3 views (last 30 days)
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

Accepted Answer

Cris LaPierre
Cris LaPierre on 26 Sep 2021
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

More Answers (1)

Paul
Paul on 26 Sep 2021
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.

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!