ODE45 diverges on specific initial conditions

hi guys,

I'm trying to run the code added below and it seems to work just fine, the only problem is that if I set my initial conditions to be [0 0 0 0] I get y to be a matrix of NaN.

When I change the initial conditions to [0.001 0 0 0] (meaning changing only the first initial condition) it works just fine. I'm guessing somwhere in the odeFun something is divided by the first initial condition.

Anybody knows if this is solvable somehow?

Thanks!

clc,clear, close

% Define symbolic variables
syms theta1(t) theta2(t) delta_t1 delta_t2

% Define parameters
m1 = 0.15; m2 = 0.15;
R = 0.2; g = 9.81;
d = 0.15; a1 = 0.05; a2 = 0.05;
k = 1.5; l0 = 0.01;
c_l = 0.1; c_0 = 0.001;
M1 = 0; M2 = 0; p = 0;

% Define expressions
theta1_dot = diff(theta1, t);
theta2_dot = diff(theta2, t);
r1 = [a1*sin(theta1), a1*cos(theta1)];
r2 = [d+a2*sin(theta2), a2*cos(theta2)];
l = norm(r2-r1);
e_l = (r2-r1)/l;

T = 0.5*m1*(R*theta1_dot)^2 + 0.5*m2*(R*theta2_dot)^2;
V = 0.5*k*(l-l0)^2 + m1*g*R*(1-cos(theta1)) + m2*g*R*(1-cos(theta2));
D = 0.5*(c_l*(diff(l, t))^2 + c_0*(theta1_dot)^2 + c_0*(theta2_dot)^2);

Fn1 = k*(l-l0)*e_l;
Fn2 = -Fn1;
drn1 = diff(r1,theta1)*delta_t1;
drn2 = diff(r2,theta2)*delta_t2;

w4 = Fn1*drn1.';
w5 = Fn2*drn2.';
w1 = M1*delta_t1;
w2 = M2*delta_t2;
w3 = 0.5*p*cos(theta1)*(R^2-a1^2)*delta_t1;
W = w1+w2+w3+w4+w5;
Q1 = diff(W,delta_t1);
Q2 = diff(W,delta_t2);
L = T-V;

% Define equations
eq1 = Q1 == diff(diff(L,theta1_dot),t)-diff(L,theta1)+diff(D,theta1_dot);
eq2 = Q2 == diff(diff(L,theta2_dot),t)-diff(L,theta2)+diff(D,theta2_dot);

% Solve ODEs
[F,~] = odeToVectorField(eq1, eq2);
odeFun = matlabFunction(F, 'Vars', {t,'Y'});
[t, y] = ode45(odeFun,[0 100],[0.0001 0 0 0]);
plot(t, y);
legend({'$\theta1$', '$\dot{\theta1}$', '$\theta2$', '$\dot{\theta2}$'},'FontSize', 16,'Interpreter', 'latex','Location', 'best');

[James Tursa]

clc,clear, close

% Define symbolic variables
syms theta1(t) theta2(t) delta_t1 delta_t2

% Define parameters
m1 = 0.15; m2 = 0.15;
R = 0.2; g = 9.81;
d = 0.15; a1 = 0.05; a2 = 0.05;
k = 1.5; l0 = 0.01;
c_l = 0.1; c_0 = 0.001;
M1 = 0; M2 = 0; p = 0;

% Define expressions
theta1_dot = diff(theta1, t);
theta2_dot = diff(theta2, t);
r1 = [a1*sin(theta1), a1*cos(theta1)];
r2 = [d+a2*sin(theta2), a2*cos(theta2)];
l = norm(r2-r1);
e_l = (r2-r1)/l;

T = 0.5*m1*(R*theta1_dot)^2 + 0.5*m2*(R*theta2_dot)^2;
V = 0.5*k*(l-l0)^2 + m1*g*R*(1-cos(theta1)) + m2*g*R*(1-cos(theta2));
D = 0.5*(c_l*(diff(l, t))^2 + c_0*(theta1_dot)^2 + c_0*(theta2_dot)^2);

Fn1 = k*(l-l0)*e_l;
Fn2 = -Fn1;
drn1 = diff(r1,theta1)*delta_t1;
drn2 = diff(r2,theta2)*delta_t2;

w4 = Fn1*drn1.';
w5 = Fn2*drn2.';
w1 = M1*delta_t1;
w2 = M2*delta_t2;
w3 = 0.5*p*cos(theta1)*(R^2-a1^2)*delta_t1;
W = w1+w2+w3+w4+w5;
Q1 = diff(W,delta_t1);
Q2 = diff(W,delta_t2);
L = T-V;

% Define equations
eq1 = Q1 == diff(diff(L,theta1_dot),t)-diff(L,theta1)+diff(D,theta1_dot);
eq2 = Q2 == diff(diff(L,theta2_dot),t)-diff(L,theta2)+diff(D,theta2_dot);

% Solve ODEs
[F,~] = odeToVectorField(eq1, eq2)

F = 

odeFun = matlabFunction(F, 'Vars', {t,'Y'})

odeFun = function_handle with value:
    @(t,Y)[Y(2);sin(Y(1)).*(-9.81e+2./2.0e+1)-Y(2)./6.0+((abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*1.0./sqrt((sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)).*((cos(conj(Y(1))).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1))./2.0e+1+(cos(Y(1)).*(sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1))./2.0e+1)-abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*((sin(conj(Y(1))).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))./2.0e+1+(sin(Y(1)).*(cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1))./2.0e+1).*1.0./sqrt((cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))).*(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*1.0./sqrt((cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1)).*((cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*((conj(Y(2)).*sin(conj(Y(1))))./2.0e+1-(conj(Y(4)).*sin(conj(Y(3))))./2.0e+1)+(cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*((sin(Y(1)).*Y(2))./2.0e+1-(sin(Y(3)).*Y(4))./2.0e+1))-abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*1.0./sqrt((sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)).*(((conj(Y(2)).*cos(conj(Y(1))))./2.0e+1-(conj(Y(4)).*cos(conj(Y(3))))./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)+((cos(Y(1)).*Y(2))./2.0e+1-(cos(Y(3)).*Y(4))./2.0e+1).*(sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1))).*(2.5e+1./6.0))./(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2)-(abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*1.0./sqrt((sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)).*((cos(conj(Y(1))).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1))./2.0e+1+(cos(Y(1)).*(sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1))./2.0e+1)-abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*((sin(conj(Y(1))).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))./2.0e+1+(sin(Y(1)).*(cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1))./2.0e+1).*1.0./sqrt((cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))).*(sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2)-1.0./1.0e+2).*1.0./sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*1.25e+2+sin(Y(1)).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*1.0./sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(3.0./2.0)-3.0./2.0e+2).*(2.5e+1./3.0)-cos(Y(1)).*1.0./sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(3.0./2.0)-3.0./2.0e+2).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*(2.5e+1./3.0);Y(4);sin(Y(3)).*(-9.81e+2./2.0e+1)-Y(4)./6.0-((abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*1.0./sqrt((sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)).*((cos(conj(Y(3))).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1))./2.0e+1+(cos(Y(3)).*(sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1))./2.0e+1)-abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*((sin(conj(Y(3))).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))./2.0e+1+(sin(Y(3)).*(cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1))./2.0e+1).*1.0./sqrt((cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))).*(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*1.0./sqrt((cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1)).*((cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*((conj(Y(2)).*sin(conj(Y(1))))./2.0e+1-(conj(Y(4)).*sin(conj(Y(3))))./2.0e+1)+(cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*((sin(Y(1)).*Y(2))./2.0e+1-(sin(Y(3)).*Y(4))./2.0e+1))-abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*1.0./sqrt((sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)).*(((conj(Y(2)).*cos(conj(Y(1))))./2.0e+1-(conj(Y(4)).*cos(conj(Y(3))))./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)+((cos(Y(1)).*Y(2))./2.0e+1-(cos(Y(3)).*Y(4))./2.0e+1).*(sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1))).*(2.5e+1./6.0))./(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2)+(abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*1.0./sqrt((sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1)).*((cos(conj(Y(3))).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1))./2.0e+1+(cos(Y(3)).*(sin(conj(Y(1)))./2.0e+1-sin(conj(Y(3)))./2.0e+1+3.0./2.0e+1))./2.0e+1)-abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*((sin(conj(Y(3))).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))./2.0e+1+(sin(Y(3)).*(cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1))./2.0e+1).*1.0./sqrt((cos(conj(Y(1)))./2.0e+1-cos(conj(Y(3)))./2.0e+1).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1))).*(sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2)-1.0./1.0e+2).*1.0./sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*1.25e+2-sin(Y(3)).*(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).*1.0./sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(3.0./2.0)-3.0./2.0e+2).*(2.5e+1./3.0)+cos(Y(3)).*1.0./sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(sqrt(abs(cos(Y(1))./2.0e+1-cos(Y(3))./2.0e+1).^2+abs(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).^2).*(3.0./2.0)-3.0./2.0e+2).*(sin(Y(1))./2.0e+1-sin(Y(3))./2.0e+1+3.0./2.0e+1).*(2.5e+1./3.0)]

%[t, y] = ode45(odeFun,[0 100],[0.0001 0 0 0]);
%plot(t, y);
%legend({'$\theta1$', '$\dot{\theta1}$', '$\theta2$', '$\dot{\theta2}$'},'FontSize', 16,'Interpreter', 'latex','Location', 'best');

Running this without the ode45() call and just examining F, you can see that there will be divide-by-0 in some of the terms, e.g., sigma7 and sigma10.

[roy]

is there anyway to make it work?

when I calculate the system analytcally theres no such problem.

[Torsten]

How do you "calculate the system analytically" ?

[roy]

pretty much the same way only manually on paper.

based on the following system:

[Torsten]

You want to tell us that you solve the above differential equations resulting in eq1 and eq2 analytically with paper and pencil ?

[roy]

hahaha no, sorry for the confusion.

I only got the differential equations resulting in eq1 and eq2 analytically.

after that I put the very long diffrential equations in Matlab and it worked for every initial condition.

I just wanted to get to those equations using matlab instead of writing all of them manually.

[Torsten]

If it works for your deduced equations and it doesn't work for eq1 and eq2, the two must be different. Can you directly compare your equations and eq1 and eq2 ? At least by evaluating them for certain arguments ?

[Sam Chak]

@roy, could you kindly demonstrate the analytical derivation of the Lagrange equations (Eq.1 and Eq.2)? This will allow us to compare them with the results obtained from odeToVectorField().