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 
