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');
댓글 수: 8
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');
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
2024년 2월 28일
Torsten
2024년 2월 28일
How do you "calculate the system analytically" ?
roy
2024년 2월 28일
pretty much the same way only manually on paper.
based on the following system:
Torsten
2024년 2월 28일
You want to tell us that you solve the above differential equations resulting in eq1 and eq2 analytically with paper and pencil ?
roy
2024년 2월 28일
답변 (0개)
카테고리
도움말 센터 및 File Exchange에서 Programming에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
