I get this error:"This system does not seem to be linear."

6 visualizzazioni (ultimi 30 giorni)
When I try this code for simple pendulum, I get results:
syms pt(t) th(t)
m=1; l=0.5; g=9.81;
e1= diff(th)*(m*l^2)==pt;
e2= diff(pt)==-m*g*l*sin(th);
vars = [pt(t); th(t)];
V = odeToVectorField([e1,e2]);
M = matlabFunction(V, 'vars', {'t','Y'});
interval = [0 5];
y0 = [0; pi/4];
ySol = ode45(M,interval,y0);
tValues = linspace(interval(1),interval(2),1000);
a= deval(ySol,tValues,1)/(m*l^2);
plot(tValues,a)
But when I use it for triple pendulum, it gives error. Couldn't solve it. Sorry if it's simple to figure out. Really new here.
syms theta1(t) theta2(t) theta3(t) p1(t) p2(t) p3(t)
m1=1; m2=1; m3=1; l1=1; l2=1;
l3=1; g=9.81; tau1=0; tau2=0; tau3=0;
I1=0; I2=0; I3=0;
e1= diff(theta1)*(I1+(m1+m2+m3)*l1^2)...
+diff(theta2)*(m2+m3)*l1*l2*cos(theta1-theta2)+...
diff(theta3)*m3*l1*l3*cos(theta1-theta3)==p1;
e2= diff(theta1)*(m2+m3)*l1*l2*cos(theta1-theta2)+...
diff(theta2)*(I2+(m2+m3)*l2^2)+...
diff(theta3)*m3*l2*l3*cos(theta2-theta3)==p2;
e3= diff(theta1)*m3*l1*l3*cos(theta1-theta3)+...
diff(theta2)*m3*l2*l3*cos(theta2-theta3)+...
diff(theta3)*(I3+m3*l3^2)==p3;
e4= diff(p1)== tau1-tau2-(m2+m3)*diff(theta1)*diff(theta2)*sin(theta1-theta2)...
-m3*diff(theta1)*diff(theta3)*l1*l3*sin(theta1-theta3)...
-(m1+m2+m3)*g*l1*cos(theta1);
e5= diff(p2)==tau2-tau3+(m2+m3)*diff(theta1)*theta2*l1*l2*sin(theta1-theta2)...
-m3*diff(theta2)*diff(theta3)*l2*l3*sin(theta2-theta3)...
-(m2+m3)*g*l2*cos(theta2);
e6= diff(p3)==tau3+(m3)*diff(theta1)*diff(theta3)*l1*l3*sin(theta1-theta3)...
+m3*diff(theta2)*diff(theta3)*l2*l3*sin(theta2-theta3)...
-m3*g*l3*cos(theta3);
vars= [theta1(t);theta2(t);theta3(t);p1(t);p2(t);p3(t)];
V = odeToVectorField([e1,e2,e3,e4,e5,e6]);
M = matlabFunction(V,'vars', {'t','Y'});
I have get that error in simple pendulum too, it was pt==diff(th)*(m*l^2), then I put the pt to the end, and it's solved. In triple pendulum I tried leaving diff(theta1) alone didn't work, tried to this code too, but nothing changed. Original equations are:
Adsız1.png
Adsız2.png

Risposta accettata

Torsten
Torsten il 4 Gen 2019
The product of differentials in your equations (diff(theta1)*diff(theta3), e.g.) makes it impossible to use ODE45.
I don't know if it can be applied directly, but ODE15I is the correct solver to use in this case.
Best wishes
Torsten.
  4 Commenti
Torsten
Torsten il 4 Gen 2019
function main
y0 = [pi/2; pi/2; pi/2; 0; 0; 0];
yp0=[0; 0; 0; 0; 0; 0;];
[y0,yp0] = decic(@odennotfunatall,0,y0,[pi/2 pi/2 pi/2 0 0 0],yp0,[]);
[t,y] = ode15i(@odennotfunatall,[0 5],y0,yp0);
plot(t,y)
end
function hell2 = odennotfunatall(~,y,yp)
m1=1; m2=1; m3=1; l1=1; l2=1;
l3=1; g=9.81; tau1=0; tau2=0; tau3=0; I1=0; I2=0; I3=0;
hell2=zeros(6,1);
hell2(1)=yp(1)*(I1+(m1+m2+m3)*l1^2)...
+yp(2)*(m2+m3)*l1*l2*cos(y(1)-y(2))+...
yp(3)*m3*l1*l3*cos(y(1)-y(3))-y(4);
hell2(2)=yp(1)*(m2+m3)*l1*l2*cos(y(1)-y(2))+...
yp(2)*(I2+(m2+m3)*l2^2)+...
yp(3)*m3*l2*l3*cos(y(2)-y(3))-y(5);
hell2(3)=yp(1)*m3*l1*l3*cos(y(1)-y(3))+...
yp(2)*m3*l2*l3*cos(y(2)-y(3))+...
yp(3)*(I3+m3*l3^2)-y(6);
hell2(4)=-yp(4)+tau1-tau2-(m2+m3)*yp(1)*yp(2)*sin(y(1)-y(2))...
-m3*yp(1)*yp(3)*l1*l3*sin(y(1)-y(3))...
-(m1+m2+m3)*g*l1*cos(y(1));
hell2(5)=-yp(5)+tau2-tau3+(m2+m3)*yp(1)*yp(2)*l1*l2*sin(y(1)-y(2))...
-m3*yp(2)*yp(3)*l2*l3*sin(y(2)-y(3))...
-(m2+m3)*g*l2*cos(y(2));
hell2(6)=-yp(6)+tau3+(m3)*yp(1)*yp(3)*l1*l3*sin(y(1)-y(3))...
+m3*yp(2)*yp(3)*l2*l3*sin(y(2)-y(3))...
-m3*g*l3*cos(y(3));
end
Arda Nova
Arda Nova il 4 Gen 2019
Yeah the zeros, couldn't find a way to put it, thank you very much.

Accedi per commentare.

Più risposte (1)

madhan ravi
madhan ravi il 3 Gen 2019
Just follow the same way showed in your previous question?
  3 Commenti
madhan ravi
madhan ravi il 3 Gen 2019
No problem , please recheck your equations take your time there are only four equations whereas you have 6 equations in your code.
Arda Nova
Arda Nova il 3 Gen 2019
Oh sorry, upper picture is streched out. At the right of "(60)" (on the picture), there is two more equations.

Accedi per commentare.

Categorie

Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by