Solving for the Motion of a Spring Pendulum

[hutchphd]

My apologies but I have lost track of what you are trying to do here. Where is r on this drawing? Have you included the parallel axis theorem in your consideration? I'm lost. Perhaps someone else understands.

 

[erobz]

My apologies but I have lost track of what you are trying to do here. Where is r on this drawing? Have you included the parallel axis theorem in your consideration? I'm lost. Perhaps someone else understands.

That last mechanism was just made up to illustrate my concern. $r$ would be the distance from the pivot to the center of the "wheel" on that drawing. And the $I$ for this system would use the parallel axis theorem, its just not clear how $I$ factors into the terms for $\tau$ above.

See, when I first was shown:

$$ L = m r^2 \dot \theta $$

the term $2 r \dot r$ came from differentiating $L$, and for a point mass the term $2 r \dot r = \dot I$, but I don't see how that is the case for the actual $I$ which might be any function.

The only thing is I would think it applies to each mass in the system individually, so in that sense the definition $\dot L = \frac{d}{dt} \left( m r^2 \dot \theta \right) \neq \frac{d}{dt}\left( I \dot \theta \right)$

I'm sorry if I am bad at explaining myself, but I am trying. It also may be the case that I have talked myself out of rationality. I'll probably wake up and see what mistake I am making. Thanks for your help.

 

[erobz]

I think this is the crux of my issue. It isn't apparent to me that the following definitions are equivalent:

$$\tau = \frac{d}{dt} \left( \mathbf{r} \times m \mathbf{v}\right) $$

$$ \tau = \frac{d}{dt} \left( I \dot \theta \right) \mathbf{u}$$

Again, $\mathbf{u}$ is just a unit vector placeholder for whatever coordinates you want to work in ( I've been scolded in PF for equating scalars and vectors previously! )

If I'm wrong, I'd like to see how?

 

[erobz]

Just to get back on track:

$(2)$ already appears to agree with the Wiki:

$$ - m \ddot r +mg \cos \theta -k\left( r - l_o \right) + mr {\dot \theta}^2 = 0 \tag{2}$$

The only thing is there is a sign difference in between the external torque terms between $(1)$ and the Wiki derivation if I take $\mathbf{r} \times m \mathbf{v} = mr^2 \dot \theta$. I'm thinking that has to come from the cross product based on my assumed convention. So $(1)$ should actually be:

$$ \sum \tau = \frac{d}{dt}\left( -m r^2 \dot \theta \right)$$

$$ rmg \sin \theta = -m \left( 2 r \dot r \dot \theta + r^2 \ddot \theta \right)$$

$$ r^2 \ddot \theta + 2 r \dot r \dot \theta + rg \sin \theta = 0 \tag{1} $$

https://en.wikipedia.org/wiki/Elastic_pendulum#Lagrangian

 

[bob012345]

I think this is the crux of my issue. It isn't apparent to me that the following definitions are equivalent:

$$\tau = \frac{d}{dt} \left( \mathbf{r} \times m \mathbf{v}\right) $$

$$ \tau = \frac{d}{dt} \left( I \dot \theta \right) \mathbf{u}$$

Again, $\mathbf{u}$ is just a unit vector placeholder for whatever coordinates you want to work in ( I've been scolded in PF for equating scalars and vectors previously! )

If I'm wrong, I'd like to see how?

I think it would be helpful to go back to fundamental definitions. Torque is defined as
$$ \large\vec{\tau} = \mathbf{r} \times \mathbf{F} $$

Expanding your first equation above
$$\large\vec{\tau} = \frac{d}{dt} \left( \mathbf{r} \times m \mathbf{v}\right) $$
$$\large\vec{\tau} = \frac{d}{dt} \left( \mathbf{r} \times \mathbf{p}\right) $$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{r} \times \mathbf{p} + \mathbf{r} \times \frac{d}{dt} \mathbf{p}$$

$$\large\vec{\tau} = \mathbf{v} \times \mathbf{p} + \mathbf{r} \times \mathbf{F}$$
but the first term is zero since $\mathbf{v}$ and $\mathbf{p}$ are parallel leaving

$$ \large\vec{\tau} = \mathbf{r} \times \mathbf{F} $$

Your second equation

$$\tau = \frac{d}{dt} \left( I \dot \theta \right) \mathbf{u}$$

is problematic if the unit vector $\mathbf{u}$ is in the direction of $\mathbf{v}$. It would be better to start with the definition of $\mathbf{L} = \mathbf{r} \times \mathbf{p}$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{L} = \frac{d}{dt} (I \mathbf{\vec{\omega}})$$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{L} = I\frac{d}{dt}\mathbf{\vec{\omega}} + \mathbf{\vec{\omega}}\frac{d}{dt}(I) $$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{L} = I\mathbf{\vec{\alpha}} + \mathbf{\vec{\omega}}\frac{d}{dt}(I ) $$
where $\mathbf{\omega}$ is a vector and $I$ a scaler.

Both ways are equivalent. The second term is not necessarily zero if $I$ changes as it does in this problem.

https://en.wikipedia.org/wiki/Torque

 

[erobz]

I think it would be helpful to go back to fundamental definitions. Torque is defined as
$$ \large\vec{\tau} = \mathbf{r} \times \mathbf{F} $$

Expanding your first equation above
$$\large\vec{\tau} = \frac{d}{dt} \left( \mathbf{r} \times m \mathbf{v}\right) $$
$$\large\vec{\tau} = \frac{d}{dt} \left( \mathbf{r} \times \mathbf{p}\right) $$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{r} \times \mathbf{p} + \mathbf{r} \times \frac{d}{dt} \mathbf{p}$$

$$\large\vec{\tau} = \mathbf{v} \times \mathbf{p} + \mathbf{r} \times \mathbf{F}$$
but the first term is zero since $\mathbf{v}$ and $\mathbf{p}$ are parallel leaving

$$ \large\vec{\tau} = \mathbf{r} \times \mathbf{F} $$

Your second equation

$$\tau = \frac{d}{dt} \left( I \dot \theta \right) \mathbf{u}$$

is problematic if the unit vector $\mathbf{u}$ is in the direction of $\mathbf{v}$. It would be better to start with the definition of $\mathbf{L} = \mathbf{r} \times \mathbf{p}$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{L} = \frac{d}{dt} (I \mathbf{\vec{\omega}})$$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{L} = I\frac{d}{dt}\mathbf{\vec{\omega}} + \mathbf{\vec{\omega}}\frac{d}{dt}(I) $$

$$\large\vec{\tau} = \frac{d}{dt} \mathbf{L} = I\mathbf{\vec{\alpha}} + \mathbf{\vec{\omega}}\frac{d}{dt}(I ) $$
where $\mathbf{\omega}$ is a vector and $I$ a scaler.

Both ways are equivalent. The second term is not necessarily zero if $I$ changes as it does in this problem.

https://en.wikipedia.org/wiki/Torque

The thing is we have to derive $I = \int r^2 dm$ into the definition:

$$ \vec{ \tau} = \vec{r} \times \vec{F}$$

With all due respect: I believe the math, and I believe that I am wrong about this issue, but the $I$ just seems to have appeared in your explanation and the Wiki article for that matter.

 

[bob012345]

The thing is we have to derive $I = \int r^2 dm$ into the definition:

$$ \vec{ \tau} = \vec{r} \times \vec{F}$$

With all due respect: I believe the math, and I believe that I am wrong about this issue, but the $I$ just seems to have appeared in your explanation and the Wiki article for that matter.

It comes from the definition $\vec{L} = I \vec{\omega}$.

 

[erobz]

I think I found what I was looking for in my Dynamics textbook:

Thanks for keeping me honest! Sorry if I gave people a headache.

 

[bob012345]

Thanks for keeping me honest! Sorry if I gave people a headache.

Are you still planning on solving this problem?

 

[erobz]

Are you still planning on solving this problem?

I've probably tortured everyone here enough!

 

[bob012345]

I've probably tortured everyone here enough!

Nonsense. I think people here enjoy helping others work through problems. If you enjoy this problem and want to see it through then don't worry. I think you can get the equations now or rationalize how they were derived in the link @hutchphd gave earlier for the elastic pendulum. I know the equations boil down to two linked ODE's but I am not sure if they can be separated completely by some variable transform.

 

[erobz]

Just to get back on track:

$(2)$ already appears to agree with the Wiki:

$$ - m \ddot r +mg \cos \theta -k\left( r - l_o \right) + mr {\dot \theta}^2 = 0 \tag{2}$$

The only thing is there is a sign difference in between the external torque terms between $(1)$ and the Wiki derivation if I take $\mathbf{r} \times m \mathbf{v} = mr^2 \dot \theta$. I'm thinking that has to come from the cross product based on my assumed convention. So $(1)$ should actually be:

$$ \sum \tau = \frac{d}{dt}\left( -m r^2 \dot \theta \right)$$

$$ rmg \sin \theta = -m \left( 2 r \dot r \dot \theta + r^2 \ddot \theta \right)$$

$$ r^2 \ddot \theta + 2 r \dot r \dot \theta + rg \sin \theta = 0 \tag{1} $$

https://en.wikipedia.org/wiki/Elastic_pendulum#Lagrangian

I think I have the proper equations now, but I don't know how to formally attack this problem.

 

[bob012345]

I think I have the proper equations now, but I don't know how to formally attack this problem.

It seems a highly non-linear set of coupled equations. I wonder is some variable transformation might simplify them? Or writing the equations in terms of x and y coordinates?

Or you could guess at the form of a general solution and solve for the constants.

 

[Baluncore]

Or you could guess at the form of a general solution and solve for the constants.

I have the recurring nightmare, of a Phase Locked Loop that fails to lock, but bounces around forever in a noisy, quasi-chaotic pattern. The probability of such a PLL is near zero, but given two lightly coupled PLLs, it becomes much more likely. That is what you have here.

 

[erobz]

It seems a highly non-linear set of coupled equations. I wonder is some variable transformation might simplify them? Or writing the equations in terms of x and y coordinates?

This isn't saying much, but that substitution is probably something I couldn't find in my wildest dreams!

Or you could guess at the form of a general solution and solve for the constants.

I think one needs to have much more experience than myself to do this. Sorry for my defeatist attitude!

 

[hutchphd]

The paper in Physica D (which I can access via alumna privileges I guess) is a pretty complete numerical treatment. The abstract is is descriptive worth a read) and echoes this behavior. Paper is nicely filled with trajectories. Rich but not simple.

 

[erobz]

The paper in Physica D (which I can access via alumna privileges I guess) is a pretty complete numerical treatment. The abstract is is descriptive worth a read) and echoes this behavior. Paper is nicely filled with trajectories. Rich but not simple.

Yeah, that abstract just dashed my hopes. It saying that parts of this are still not fully understood. I'm expecting that some aspects of it are doctoral level research in Classical Mechanics?!?

I'm letting this one slide (the floors in my house need swept and steamed)!

 

[bob012345]

The paper in Physica D (which I can access via alumna privileges I guess) is a pretty complete numerical treatment. The abstract is is descriptive worth a read) and echoes this behavior. Paper is nicely filled with trajectories. Rich but not simple.

I was able to reproduce the pattern shown in the Wikipedia page on elastic pendulums with an app that uses two generating circles but I suppose it is no surprise that the elastic pendulum has some conditions that will make orderly closed Lissajous figures and others that transition to chaotic motion whereas simple circle generators do not. This led me to think that there might be a closed form solution to this problem but now I think not.

https://www.geogebra.org/m/vRA9prdc

 

[bob012345]

Yeah, that abstract just dashed my hopes. It saying that parts of this are still not fully understood. I'm expecting that some aspects of it are doctoral level research in Classical Mechanics?!?

I'm letting this one slide (the floors in my house need swept and steamed)!

I am in a science interest group at my local Makerspace and we have purchased an analog computer called The Analog Thing to play with. I know analog computers can do double pendulums. At home I use an LTSpice simulator to run an analog simulation when I want to do analog. https://the-analog-thing.org/

 

[hutchphd]

I am in a science interest group at my local Makerspace and we have purchased an analog computer called The Analog Thing to play with.

That sounds like bunches of fun. Can you program (or interest someone else in programming) the planar elastic pendulum? How does it I/O ?

 

[bob012345]

That sounds like bunches of fun. Can you program (or interest someone else in programming) the planar elastic pendulum? How does it I/O ?

It's coming from Germany and should be here any day now . The Analog Thing (THAT) is limited so probably only does simpler systems but I checked the online manual and it seems to do several simpler chaotic equations

https://the-analog-thing.org/docs/dirhtml/rst/applications/elegant_chaos/alpaca_15/

The double pendulum requires ten integrators, 15 summers, 16 multipliers, and 17 coefficient potentiometers. Each THAT has 5 integrators, 4 summers, 2 multipliers and 8 coefficient pots which would require 8 THAT's so not doable in one but it can be done in an LTSpice simulator which of course can have an arbitrary number of elements and is free. Here is the setup

https://the-analog-thing.org/docs/dirhtml/rst/applications/double_pendulum/alpaca_21/

The elastic pendulum should be doable in LTSpice but I am a bit worried about the variables in the denominator for these equations. It may be a simple fix to just multiply by the denominator and rearrange the equations into products rather than divisions. I'd be willing to try it.

As far as outputs, it can be hooked to a scope or a PC with a simulated scope;
https://the-analog-thing.org/docs/dirhtml/rst/computing_elements/THAT_elements/outputs/

EDIT: Of course one could make more integrators, summers and inverters with some op amps and breadboards ect. to add to the system as well as buy some analog multiplier chips.

 

[bob012345]

I simulated this elastic pendulum using LTSpice. I set the problem in terms of state variables according to this paper which suggests Spice as a powerful general purpose tool for almost any dynamical system.

https://www.eng.auburn.edu/~wilambm/pap/2011/02_ijee2410ns.pdf

It's a powerful tool. I've even done quantum simulations with it. The entire code to run this system is here with V(1)=$\theta$, V(2)=$\frac{d\theta}{dt}$, V(3)=$r$, V(4)=$\frac{dr}{dt}$;

.param l=1.0, k=10, m=1 , g=10
C1 0 1 1
C2 0 2 1
C3 0 3 1
C4 0 4 1
R1 0 1 100G
R2 0 2 100G
R3 0 3 100G
R4 0 4 100G
G1 0 1 Value={V(2)}
G2 0 2 Value={-2*V(2)*V(4)/(l+V(3)) -g*sin(V(1))/(l+V(3))}
G3 0 3 Value={V(4)}
G4 0 4 Value={(l + V(3))*(V(2))**2 - (k/m)*V(3) +g*cos(V(1))}
.IC V(1)=1 V(2)=0 V(3)=0.57 V(4)=0
.tran 0.01ms 8.5s 0 0.1ms UIC
.probe
.end

The output is as close as I could get it to match the Wikipedia page (since I don't know the exact conditions I had to try and match the curves). It's not an absolute match but close except for the artifact at the start (I removed that thanks to @Baluncore ).

The Wiki page has this simulation;

Here is another one with a few changes. Looks like a pretzel;

The next step would be to do it as a circuit.

 

[Baluncore]

The next step would be to do it as a circuit.

There is not much circuit there, only five nodes.

The order of lines or blocks of text does not matter in spice, unless there is a .end, which stops the process, so depending on the order of processing, some lines may be missed. LTspice does not need a .end directive. Remove the .end and .probe directives from the block of code. Then you can start breaking up the text block, and making the circuit.

In the .tran directive, introduce a 1 ns "time to start saving data" delay. That will prevent display of the initial diagonal step transient.

 

[bob012345]

There is not much circuit there, only five nodes.

The order of lines or blocks of text does not matter in spice, unless there is a .end, which stops the process, so depending on the order of processing, some lines may be missed. LTspice does not need a .end directive. Remove the .end and .probe directives from the block of code. Then you can start breaking up the text block, and making the circuit.

In the .tran directive, introduce a 1 ns "time to start saving data" delay. That will prevent display of the initial diagonal step transient.

By circuit I meant using a series of integrators, summers, inverters, multipliers and such as is done in analog computers not converting this state variable format directly to a circuit diagram format with programmable conductances. Of course the ease with which one can run virtually any dynamical system this way makes it seem almost like cheating. Here is an example of a simple damped mass-spring oscillator in analog circuit form;

Thanks for the tip on the delay to remove the artifact!

 

[Baluncore]

By circuit I meant using a series of integrators, summers, inverters, multipliers and such as is done in analog computers not converting this state variable format directly to a circuits diagram format.

I guess you want to eliminate all G and B sources, along with expression evaluation.
The reason spice can be fast here, is that the number of state variables, or nodes, can be minimised to four and a common ground. This is really the simulation of a hybrid computer, the currents integrated in the virtual capacitors are generated directly from the digital equations, not by simulation of the many op-amps with feedback, needed to make multipliers, and a divider, with coefficient resistors to convert the voltage terms to currents. That will introduce a great number of nodes into the sparse matrix, which will slow it all down.

The fastest LTspice circuit simulation would use an arbitrary behavioral current source, B, in place of a G, (which has only one differential control input).
B1 0 1 I={ V(2) }
C1 0 1 1
R1 0 1 1T
B2 0 2 I={( -2*V(2)*V(4) -g*sin(V(1) ) )/(l+V(3)) }
C2 0 2 1
R2 0 2 1T
B3 0 3 I={ V(4) }
C3 0 3 1
R3 0 3 1T
B4 0 4 I={ (l + V(3))*(V(2))**2 -(k/m)*V(3) +g*cos(V(1)) }
C4 0 4 1
R4 0 4 1T

 

[TSny]

As a check, I used Mathematica to solve the differential equations. I get the same results as @bob012345 for his two examples:

The axes are scaled a little differently.

Using software for electric circuits to numerically solve the equations of motion of a mechanical system is really neat.

 

[erobz]

Now that see all these interesting technical approaches that are completely over my head, I feel the urge to go cave man on this. Maybe...

I want to try to outline what I intend to do in full, so it's more apparent if there are glaringly obvious "no way is that going to work" errors before proceeding with the programming part of it.

Begining with the equations of motion:

$$ r^2 \ddot \theta + 2 r \dot r \dot \theta + rg \sin \theta = 0 \tag{1}$$

$$ \ddot r + \left( \frac{k}{m} - {\dot \theta}^2 \right) r = g\cos \theta +\frac{k l_o}{m} \tag{2}$$

So I want to work with $(2)$ treating the coefficients as constant. This implies the following form:

$$ \ddot r + \beta^2 r = \gamma $$

Where;

$$ \beta = \sqrt{ \frac{k}{m} - {\dot \theta}^2} $$

$$ \gamma = g\cos \theta +\frac{k l_o}{m} $$

That equation has the solution:

$$ r(t) = C_1 \sin( \beta t ) + C_2 \cos ( \beta t ) + \frac{\gamma}{\beta^2} $$

Then the objective is to solve for $C_1$and $C_2$ at the beginning of each time step $t_o$ with conditions that $r(t_o) = r_o$ and $\dot r (t_o) = \dot {r_o}$. The subcript $o$ will be used to refer to the value of each parameter at the beginning of each time step. Applying those conditions gives the following system:

$$ \begin{align} A C_1 + B C_2 &= r_o - \frac{\gamma_o}{\beta_o^2} \tag*{} \\ \quad \tag*{} \\ \beta_o B C_1 - \beta_o A C_2 &= \quad \dot {r_o} \tag*{} \end{align}$$

Where:

$$ A = \sin ( \beta_o t_o ) $$

$$ B = \cos ( \beta_o t_o ) $$

After obtaining a solution for $C_1, C_2$, I can find $r, \dot r$ at $\theta_o$. Then I will use the result above as a substition into $(1)$ isolated for the angular acceleration which will be treated as a constant over the small time interval $\Delta t$:

$$ \ddot \theta_o = -\frac{1}{r_o} \left( g \sin \theta_o + 2 \dot r_o \dot \theta_o \right) $$

From that it follows that:

$$ \dot \theta ( t + \Delta t ) = \dot \theta_o + \ddot \theta_o \Delta t $$

and

$$ \theta ( t + \Delta t ) = \theta_o + \dot \theta_o \Delta t + \frac{1}{2} \ddot \theta_o { \Delta t }^2 $$

That is my full plan, could it work?

 

[bob012345]

I want to try to outline what I intend to do in full, so it's more apparent if there are glaringly obvious "no way is that going to work" errors before proceeding with the programming part of it.

...

That is my full plan, could it work?

You seem to want to invent a new way to numerically solve coupled ODE's. You can try and see if it works but why not spend a little time first investigating basic numerical algorithms and see what is available?

 

[erobz]

You seem to want to invent a new way to numerically solve coupled ODE's. You can try and see if it works but why not spend a little time first investigating basic numerical algorithms and see what is available?

Because I think that stuff is complex. Just saying the phrase "Solution Techniques for Systems of Coupled Non-Linear ODE's" gives me chills! I would have to take multiple courses in Differential Equations and Linear Algebra before I begin to understand it. It's like a surgeon handing a toddler a scalpel and saying, "remove the malignant tissue attached to my patients' heart". In theory they could do it!

For people that are really smart (like you guys), It seems to me it's much simpler for them to look at the math and say, "here is where your first bone headed mistake was"... Something like that has already happened in this thread multiple times!

P.S. I only ask because I'm just about as good at programming as I am at solving Physics problems. So it might take me a significant amount of time to realize its crap.

 

[bob012345]

Because I think that stuff is complex. Just saying the phrase "Solution Techniques for Systems of Coupled Non-Linear ODE's" gives me chills! I would have to take multiple courses in Differential Equations and Linear Algebra before I begin to understand it. It's like a surgeon handing a toddler a scalpel and saying, "remove the malignant tissue attached to my patients' heart". I theory they could do it!

It is complicated but you don't need to get a doctorate first to do this problem. What you need perhaps is advice on what techniques to start with and where to go to learn them quickly.

For people that are really smart (like you guys), It seems to me it's much simpler for them to look at the math and say, "here is where your first bone headed mistake was"... Something like that has already happened in this thread multiple times!

You yourself derived the equations of motion above. I doubt any random stranger off the street could do that. Now you want help numerically solving them. We all have different backgrounds and levels of experience but yes, the mentors and advisors here are very smart and good at what they do. Some members like you and I are still learning.

P.S. I only ask because I'm just about as good at programming as I am at solving Physics problems. So it might take me a significant amount of time to realize its crap.

If your bottom line is to solve this problem numerically then it is well worth the trouble to learn a little of what has been proven to be useful and valuable.