ODE45 producing NaN values for a set of differential equations

Question

Mishal Mohanlal 2021-9-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1461094-ode45-producing-nan-values-for-a-set-of-differential-equations

评论： Mishal Mohanlal 2021-9-28

采纳的回答： Cris LaPierre

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Cris LaPierre 2021-9-26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1461094-ode45-producing-nan-values-for-a-set-of-differential-equations#answer_795599

在 MATLAB Online 中打开

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