Solving system of second order differential equations with ode45

4 vues (au cours des 30 derniers jours)
Thomas Holmes
Thomas Holmes le 28 Avr 2019
Commenté : Star Strider le 29 Avr 2019
I have two second equations
theta1''= [-29.4sin(theta1)- 9.8sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(theta1'^(2)*cos(theta1-theta2)]/(3-cos(2*theta1-2*theta2))
theta2''=[2*sin(theta1-theta2)[(2*theta1')+19.6cos(theta1)+theta2'*cos(theta1-theta2)]]/(3-cos(2*theta1-2*theta1))
I think these should be written as a system of 4 first order equations, recast as a matrix and put into ode45 but I cannot figure out hwo to write these equatuons as 4 first first order due to the trig functions. I think I know to use ode45 once I have them in this form. Any help on how to rewrite these would be appreaciated.
  2 commentaires
Thomas Holmes
Thomas Holmes le 29 Avr 2019
These theta1 and theta2 values represent the angles between the vertical and two point masses on a double pendulum, how can I use theta1 and theta2 to plot the path the bottom pendulum takes?
Star Strider
Star Strider le 29 Avr 2019
The angles will be the ‘y’ output of your ode45 call:
[t,y] = ode45(odesfcn, tspan, theta0);
With those and the geometry of your system, you can plot the trajectories of the point masses. You would have to model the masses as functions of the angles.
That is the best I can do in terms of a description.

Connectez-vous pour commenter.

Réponse acceptée

Star Strider
Star Strider le 28 Avr 2019
I prefer to let the Symbolic MAth Toolbox do the ‘heavy lifting’:
syms theta1(t) theta2(t) t Y
% theta1'' = -29.4sin(theta1)- 9.8sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(theta1'^(2)*cos(theta1-theta2)]/(3-cos(2*theta1-2*theta2))
% theta2'' = 2*sin(theta1-theta2)[(2*theta1')+19.6cos(theta1)+theta2'*cos(theta1-theta2)]]/(3-cos(2*theta1-2*theta1))
Dt11 = diff(theta1);
Dt12 = diff(Dt11);
Dt21 = diff(theta2);
Dt22 = diff(Dt21);
Eqn1 = Dt12 == (-29.4*sin(theta1)- 9.8*sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(Dt11^(2)*cos(theta1-theta2)))/(3-cos(2*theta1-2*theta2));
Eqn2 = Dt22 == (2*sin(theta1-theta2)*((2*theta1')+19.6*cos(theta1)+Dt21*cos(theta1-theta2)))/(3-cos(2*theta1-2*theta1));
Eqns = [Eqn1; Eqn2];
[VF,Sbs] = odeToVectorField(Eqns);
Sbs
odesfcn = matlabFunction(VF, 'Vars',{t, Y})
producing:
Sbs =
theta2
Dtheta2
theta1
Dtheta1
odesfcn =
function_handle with value:
@(t,Y)[Y(2);-sin(Y(1)-Y(3)).*(conj(Y(3)).*2.0+cos(Y(3)).*(9.8e+1./5.0)+cos(Y(1)-Y(3)).*Y(2));Y(4);-(sin(Y(1).*2.0-Y(3)).*(4.9e+1./5.0)-sin(Y(3)).*(1.47e+2./5.0)+cos(Y(1)-Y(3)).*sin(Y(1)-Y(3)).*Y(4).^2.*2.0)./(cos(Y(1).*2.0-Y(3).*2.0)-3.0)]
or:
odesfcn = @(t,Y)[Y(2);-sin(Y(1)-Y(3)).*(conj(Y(3)).*2.0+cos(Y(3)).*(9.8e+1./5.0)+cos(Y(1)-Y(3)).*Y(2));Y(4);-(sin(Y(1).*2.0-Y(3)).*(4.9e+1./5.0)-sin(Y(3)).*(1.47e+2./5.0)+cos(Y(1)-Y(3)).*sin(Y(1)-Y(3)).*Y(4).^2.*2.0)./(cos(Y(1).*2.0-Y(3).*2.0)-3.0)];
Use the ‘Sbs’ vector in your legend call.
NOTE — I did not test this with an actual ode45 call. You may need to experiment with it to get it to work correctly.
  4 commentaires
Thomas Holmes
Thomas Holmes le 28 Avr 2019
Thank you!
Star Strider
Star Strider le 28 Avr 2019
As always, my pleasure!

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by