A Lagrangian for a double spring pendulum connected through a rigid bar

AI Thread Summary
The discussion revolves around solving a Lagrangian problem for a double spring pendulum connected by a rigid bar, with the user seeking assistance due to unexpected results in their MATLAB simulations. The primary issue noted is that the "theta2" variable consistently yields a zero solution, and the user suspects that the holonomic constraint may be the source of the problem. Suggestions include simplifying the system by reducing degrees of freedom and re-expressing the variables to achieve cleaner equations. The user has managed to get the system to solve in MATLAB but has yet to verify the results' accuracy. Overall, the conversation focuses on troubleshooting the Lagrangian formulation and improving the simulation approach.
Mishal0488
Messages
20
Reaction score
1
Hi Guys

I am looking for some guidance with regards to a Lagrangian problem I am trying to solve.
Please refer to the attached documents.

Please neglect all the forcing functions for the time being. I am currently just trying to simulate the problem using initial conditions only

I have written Matlab code to simulate the problem, however my results are not as expected. the "theta2" variable always has a 0 solution, in some cases imaginary numbers form part of the solution. using initial conditions on the theta2 variable generally causes no solution as the ODE solvers on Matlab develop tolerance errors.

Note that the holonomic constraint whereby (b+r) is the subject of the formula is developed on Matlab using the symbolic package.

I believe the source of the issues is due to the holonomic constraint.
Can anyone provide me with some guidance on how to handle the system? and/or any mistakes which I have made. Could the system be represented in a more elegant manner?

Apologies for any funny terminology or poor representation of aspects, I am a structural engineer, Lagrangian mechanics isn't something that forms part of our undergraduate education hence this has been self taught.

Matlab:
function xp = Statespace3(t,x,r,L,Cs,Cp,Fv1,Fv2,Fh,w,w2,w3,k,m1,m2)

A = (x(5)+r);
three = A*cos(x(1))*cos(x(3));
four = A*sin(x(1))*sin(x(3));

R1 = L^2*sin(x(3))^2;
R2 = -A^2*cos(x(1))^2*sin(x(3))^2;
R3 = -A^2*cos(x(3))^2*sin(x(1))^2;
R4 = -2*A*L*cos(x(3))^2*sin(x(1));
R5 = 2*A^2*cos(x(1))*cos(x(3))*sin(x(1))*sin(x(3));
R6 = 2*A*L*cos(x(1))*cos(x(3))*sin(x(3));

one = sqrt(R1+R2+R3+R4+R5+R6);
B = L*sin(x(3))-one+three+four;

three_dot = x(6)*cos(x(1))*cos(x(3))-A*x(2)*sin(x(1))*cos(x(3))-A*x(4)*cos(x(1))*sin(x(3));
four_dot = x(6)*sin(x(1))*sin(x(3))+A*x(2)*cos(x(1))*sin(x(3))+A*x(4)*sin(x(1))*cos(x(3));

R1_dot = 2*x(4)*sin(x(3))*cos(x(3))*L^2;
R2_dot = -2*x(6)*A*cos(x(1))^2*sin(x(3))^2+2*x(2)*A^2*cos(x(1))*sin(x(1))*sin(x(3))^2-2*x(4)*A^2*cos(x(1))^2*sin(x(3))*cos(x(3));
R3_dot = -2*x(6)*A*cos(x(3))^2*sin(x(1))^2+2*A^2*x(4)*cos(x(3))*sin(x(3))*sin(x(1))^2-2*A^2*cos(x(3))^2*x(2)*sin(x(1))*cos(x(1));
R4_dot = -2*x(6)*cos(x(3))^2*sin(x(1))^2+2*x(4)*A*cos(x(3))*sin(x(3))*sin(x(1))^2-2*x(2)*A*cos(x(3))^2*sin(x(1))*cos(x(1));
R5_dot = 4*x(6)*A*cos(x(1))*cos(x(3))*sin(x(1))*sin(x(3))-2*A^2*x(2)*sin(x(1))^2*cos(x(3))*sin(x(3))-2*A^2*cos(x(1))*x(4)*sin(x(3))^2+sin(x(1))+2*A^2*cos(x(1))*x(2)*cos(x(1))^2*sin(x(3))+2*A^2*cos(x(1))*sin(x(1))*x(4)*cos(x(3))^2;
R6_dot = 2*x(6)*L*cos(x(1))*cos(x(3))*sin(x(3))-2*A*L*x(2)*sin(x(1))*cos(x(3))*sin(x(3))-2*A*L*x(4)*cos(x(1))*sin(x(3))^2+2*A*L*x(4)*cos(x(1))*cos(x(3))^2;

X_dot = (R1_dot+R2_dot+R3_dot+R4_dot+R5_dot+R6_dot);

b_dot = L*x(4)*cos(x(3))+three_dot+four_dot-X_dot/(2*(one));

%theta1 coordinate
F_t1 = A*Fh*sin(w*t)*cos(x(1));
theta1_dotdot = (F_t1-2*m1*x(2)*x(6)*A-9.81*m1*sin(x(1))*A-Cp*x(2))/(m1*A^2);

%theta2 coordinate
F_t2 = B*Fh*sin(w*t)*cos(x(3));
theta2_dotdot = (F_t2-2*m2*x(4)*b_dot*B-9.81*m2*sin(x(3))*B-Cp*x(4))/(m2*B^2);

%a coordinate
Fa = Fh*sin(w*t)*sin(x(1))+Fv1*sin(w2*t);
a_dotdot = (Fa+m1*x(2)^2*A-k*x(5)+9.81*m1*cos(x(1))-Cs*x(6))/m1;

%b coordinate
b_dotdot = (m2*x(4)^2*B+k*(B-r)*9.81*m2*cos(x(3)))/m2;
%state vector
xp = [x(2);
      theta1_dotdot;
      x(4);
      theta2_dotdot
      x(5)
      a_dotdot;];
 

Attachments

Physics news on Phys.org
For the case of single spring pendulum
L=\frac{m}{2}(\dot{r}^2+r^2\dot{\theta}^2)-mgr sin\theta -\frac{k}{2}(r-r_0)^2
\frac{\partial L}{\partial q}-\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}
for q=r,##\theta## are
mr\dot{\theta}^2 - mg \sin\theta -k (r-r_0) - m \ddot{r}
-m\dot{(r^2\dot{\theta})}- mgr cos\theta
where ##r_0## is natural length of spring.
By method of Lagrangean multiplier the equations of motions are
mr_1\dot{\theta_1}^2 - mg \sin\theta_1 -k (r_1-r_0) - m \ddot{r_1}+\lambda(t)\frac{\partial C}{\partial r_1}=0
mr_2\dot{\theta_2}^2 - mg \sin\theta_2 -k (r_2-r_0) - m \ddot{r_2}+\lambda(t)\frac{\partial C}{\partial r_2}=0
-m\dot{(r_1^2\dot{\theta_1})}- mgr_1 \cos\theta_1 +\lambda(t)\frac{\partial C}{\partial \theta_1}=0
-m\dot{(r_2^2\dot{\theta_2})}- mgr_2 \cos\theta_2 +\lambda(t)\frac{\partial C}{\partial \theta_2}=0
where
C=r_1^2+r_2^2-2r_1r_2cos(\theta_2-\theta_1)+2r_2L \sin\theta_2 -2r_1L\sin\theta_1
These four equations together with
C=0
which means the two spring pendulums are connected with a bar of length L.
should give us a solution.

I hope this would help you to check your way.
 
Last edited:
  • Like
Likes Delta2 and Mishal0488
Thank you.

With regards to the formulation of the Lagrangian we have done a similar thing with the difference being that my r1 and r2 are (a+r) and (b+r) however I think both a suitable.

I don't understand your holonomic constraint, can you please explain further?

Can the problem also be solved by simply making one of the variables of the holonomic constraint the subject of the formula (lets say r1)? We can then also take the time derivative of r1, thus using substitution we can reduce the number of DOF down to 3 instead of the current 4.
 
Connection by bar is expressed as
(r_2 \cos\theta_2-r_1 \cos\theta_1)^2+(r_2 sin\theta_2+L-r_1 \sin\theta_1)^2=L^2
This is transformable to C=0 in post #2.

I have just got to make equations and not gone to how to solve them. Numerical solution might be necessary.
 
Mishal0488 said:
Can the problem also be solved by simply making one of the variables of the holonomic constraint the subject of the formula (lets say r1)?
Another approach is to take x, y coordinates of the center of the bar and angle of bar wrt horizon as variables. We reduce parameters from four to three with the constraints built in.
 
Last edited:
@anuttarasammyak, Thank you for the assistance. I managed to get the system to solve on Matlab, I have not had a chance to check the legitimacy of the results as yet.

The only difference between your suggestion and what I have is the manner in which the spring length is represented. Your method developed cleaner equations though due to there being less variables to worry about. Still can't figure out the issue with what I initially did...
 
The rope is tied into the person (the load of 200 pounds) and the rope goes up from the person to a fixed pulley and back down to his hands. He hauls the rope to suspend himself in the air. What is the mechanical advantage of the system? The person will indeed only have to lift half of his body weight (roughly 100 pounds) because he now lessened the load by that same amount. This APPEARS to be a 2:1 because he can hold himself with half the force, but my question is: is that mechanical...
Some physics textbook writer told me that Newton's first law applies only on bodies that feel no interactions at all. He said that if a body is on rest or moves in constant velocity, there is no external force acting on it. But I have heard another form of the law that says the net force acting on a body must be zero. This means there is interactions involved after all. So which one is correct?
Let there be a person in a not yet optimally designed sled at h meters in height. Let this sled free fall but user can steer by tilting their body weight in the sled or by optimal sled shape design point it in some horizontal direction where it is wanted to go - in any horizontal direction but once picked fixed. How to calculate horizontal distance d achievable as function of height h. Thus what is f(h) = d. Put another way, imagine a helicopter rises to a height h, but then shuts off all...
Back
Top