odeToVectorField takes forever to run

hi guys,

I'm trying to run the added below code to solve 2 non-linear second-order ODEs, using the odeToVectorField function.

For some reason it takes forever to run, gets stuck running the odeToVectorField function.

Any ideas how can I make it work?

Thanks!

clc,clear,close
syms theta1(t)
syms theta2(t)
syms delta_t1
syms delta_t2
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;
r1 = [a1*sin(theta1), a1*cos(theta1)];
r2 = [d+a2*sin(theta2), a2*cos(theta2)];
e_l = (r2-r1)/(norm(r2-r1));
l = norm(r2-r1);
T = (0.5*m1*(R*diff(theta1))^2)+(0.5*m1*(R*diff(theta1))^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))^2+c_0*(diff(theta1))^2+c_0*(diff(theta2))^2);
Fn1 = k*(l-l0)*e_l;
drn1 = diff(r1,theta1)*delta_t1;
Fn2 = -Fn1;
drn2 = diff(r2,theta1)*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;
theta1_dot = diff(theta1,t);
theta2_dot = diff(theta2,t);
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);
[F,S] = odeToVectorField(eq1, eq2);
odeFun = matlabFunction(F, 'Vars', {t,'Y'});
[t, y] = ode45(odeFun,[0 10],[0 0 0 0]);
plot(t,y)

[Adrián László Szemenyei]

For me the odeToVectorField function takes 0.185936 seconds and it gives back a three dimensional problem, so I had to change the init conditions in the ode45 function.

Answer 1

[Torsten]

Are you sure you want to differentiate with respect to functions of t (like theta1, theta1_dot, theta2_dot) ?

drn1 = diff(r1,theta1)*delta_t1;
drn2 = diff(r2,theta1)*delta_t2;

theta1_dot = diff(theta1,t);
theta2_dot = diff(theta2,t);
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);

And from the pure aestetic, it should be

drn2 = diff(r2,theta2)*delta_t2;

instead of

drn2 = diff(r2,theta1)*delta_t2;

[roy]

yes, theta1 and theta2 are both functions of t.

fixed drn2, thanks for noticing!

the code is written based on the following system: