I Solving for the Motion of a Spring Pendulum

[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.

 

[erobz]

What you need perhaps is advice on what techniques to start with and where to go to learn them quickly.

Well, I guess there is really no rush. I'm going to try and program the numerical solution in VBA and see what happens. If it fails, it fails. The bottom line was to learn something, and I've already done that.

Sitting down and trying to learn something completely new always sucks the life out of me. Also, I rarely have uninterrupted free time to do so (I take care of my young children...stay at home dad isn't as carefree as it sounds!), but I'll see what I can dig up and welcome any suggestions.

 

[hutchphd]

If it fails, it fails. The bottom line was to learn something, and I've already done that.

Of course one of the problems for this system is that it may be difficult to parse the numerical issues from the intrinsic instabilities of the system. While I like your bottom line, I find it a little bit counterintuitive that you cannot learn the numerical technique from others but think it easier to figure it out yourself. There is no overarching utility in reinventing the wheel, so I suggest a combination approach.

 

[erobz]

Of course one of the problems for this system is that it may be difficult to parse the numerical issues from the intrinsic instabilities of the system. While I like your bottom line, I find it a little bit counterintuitive that you cannot learn the numerical technique from others but think it easier to figure it out yourself. There is no overarching utility in reinventing the wheel, so I suggest a combination approach.

I think there is utility in it though. Sometimes I find that if I try to develop it myself, I'm positioned to better understand what others are trying to tell me about it. Other times I just get horribly lost!

I think it's been said, "Through hardship comes enlightenment" and it resonated with me.

 

[erobz]

If anyone is interested, Iv'e done some manual iteration in PTC Mathcad. I've set the spring constant $k$ to be very large and initial angle $\theta_o$ to be very small in hopes to recover the small angle approximation fixed length pendulum:

$$ \ddot \theta + \frac{g}{l_o} \theta = 0 \implies \theta = \theta_o \cos \left( \sqrt{ \frac{g}{l_o} } t\right) $$

That has a quarter period $t$ of about:

$$ t = \sqrt{ \frac{l_o}{g} } \frac{\pi}{2} \approx 0.158 \rm{s} $$

The model has a qaurter period (with a somewhat large timestep) $\Delta t = 0.01 \rm{s}$ in between $0.14 \rm{s} - 0.15 \rm{s}$.

So it is undershooting just a bit, but it has a larger timestep than I would do in the actual numerical solution, and it would be approximating:

$$ \ddot \theta + \frac{g}{l_o} \sin \theta = 0 \neq \ddot \theta + \frac{g}{l_o} \theta $$

So, is the "error" in the proper direction?

If its not too much trouble would @bob012345 be able to compare some other sets of parameters like (and this set) this for short durations of time to see were they end up relative to each other?

 

[bob012345]

If its not too much trouble would @bob012345 be able to compare some other sets of parameters like (and this set) this for short durations of time to see were they end up relative to each other?

I'd be happy to try. Give me a few minutes to set it up...

 

[erobz]

I'd be happy to try. Give me a few minutes to set it up...

Thank You. No hurry!

 

[bob012345]

Thank You. No hurry!

Here are my numbers. I had to use k=100000 which I think is effectively infinity.

time Theta (deg) Omega
0.00e-00 9.740283e-01 -1.699498e-08
1.00e-02 9.691587e-01 -1.697514e-02
2.00e-02 9.545089e-01 -3.385784e-02
3.00e-02 9.302843e-01 -5.036643e-02
4.00e-02 8.969097e-01 -6.628610e-02
5.00e-02 8.544609e-01 -8.159551e-02
6.00e-02 8.034811e-01 -9.608209e-02
7.00e-02 7.444282e-01 -1.096150e-01
8.00e-02 6.784094e-01 -1.219674e-01
9.00e-02 6.052226e-01 -1.331803e-01
1.00e-01 5.258211e-01 -1.430802e-01
1.10e-01 4.414393e-01 -1.515175e-01
1.20e-01 3.524622e-01 -1.584575e-01
1.30e-01 2.602780e-01 -1.637955e-01
1.40e-01 1.650533e-01 -1.675184e-01
1.50e-01 6.910897e-02 -1.695483e-01

 

[erobz]

I changed my $\Delta t = 0.005 \rm{s}$ because I wasn't sure what level of precision your simulator was using.

You are the black plot $\theta_b$ , mine is the red plot $\theta$

Not a perfect match, but there seems to be some aberrations that could be blurring it. I don't know...I guess It is probably worth the full treatment.

Yours has that slight offset from $1 ^\circ$ and seeing them plotted, I found that I had double defined my all my parameters on the first step. The tangent to the graph at $t=0$ should be horizontal (I can see that the red curve is not), I'll fix that.

 

[bob012345]

I changed my $\Delta t = 0.005 \rm{s}$ because I wasn't sure what level of precision your simulator was using.

You are the black plot $\theta_b$ , mine is the red plot $\theta$

View attachment 303602

Not a perfect match, but there seems to be some aberrations that could be blurring it. I don't know...I guess It is probably worth the full treatment.

I wasn't sure about degree/radian conversion as far as precision. I used exactly 0.017 rads as starting angle. Also, the last data point is 0.07 not 0.007. My timesteps less that 0.1ms. My numbers don't change much if I use 10ms as the max timestep.

 

[erobz]

I wasn't sure about degree/radian conversion as far as precision. I used exactly 0.017 rads as starting angle. Also, the last data point is 0.07 not 0.007. My timesteps less that 0.1ms. My numbers don't change much if I use 10ms as the max timestep.

Mine was whatever was store in PTC software for 1 degree. I don't think I'm going to bother to change it. Its time consuming to manually go through all the iterations to update the table.

Fixed

Now I'm feeling much better! What do you think? Maybe I should change mine though, it looks like it will dive under yours if I do.

 

[bob012345]

Mine was whatever was store in PTC software for 1 degree. I don't think I'm going to bother to change it. Its time consuming to manually go through all the iterations to update the table.

Fixed

View attachment 303608

Now I'm feeling much better! What do you think? Maybe I should change mine though, it looks like it will dive under yours if I do.

Why don't you just try carrying it through one complete cycle and see if it tracks. I can send you data for the rest of the cycle.

 

[erobz]

I'm doing it "manually" right now, the number of computations are getting too large for the program.

I've corrected the offset.

Given that its reasonably in agreement, I probably will write the numerical solution with traditional programming. Going away for the holiday weekend though. Maybe next week?

It will be neat to see if that amount of disagreement given the potential of chaotic motion will spell disaster. Then again, when I am able to decrease my timestep perhaps they will get better.

 

[erobz]

If anyone wants to examine what produced the latest graph see the attached PDFEDIT:
I just decreased my time step to 0.001 s for up to time 0.01 s and it is moving closer to your simulation. Thats a good sign. I'll let it at that for now.

 

[erobz]

I tried to program it, I think it has major theoretical issues that are tied to the assumed form of $(2)$.

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

For the solution of $2$ to be of the form:

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

the condition is:

$$ \beta = \sqrt{ \frac{k}{m} - \dot{ \theta }^2 } > 0 \quad \rm{condition 1}$$

So, it is quickly failing for any sizeable angle ( non-real under the root, violating $\rm{condition 1}$ )

how about those standard numerical techniques? Any good information?

 

[bob012345]

I tried to program it, I think it has major theoretical issues that are tied to the assumed form of $(2)$.

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

For the solution of $2$ to be of the form:

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

the condition is:

$$ \beta = \sqrt{ \frac{k}{m} - \dot{ \theta }^2 } > 0 \quad \rm{condition 1}$$

So, it is quickly failing for any sizeable angle ( non-real under the root, violating $\rm{condition 1}$ )

how about those standard numerical techniques? Any good information?

Some basic information is here;

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

But I think you really would need to put in some effort to be able to do this problem numerically from scratch. You should start with much simpler problems or the limiting cases of this problem in which the equations reduce to simpler forms. Alternatively, you can learn to use Mathematica, Matlab or some other simulator to do it for you and explore the solution space.

 

[erobz]

Some basic information is here;

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

But I think you really would need to put in some effort to be able to do this problem numerically from scratch. You should start with much simpler problems or the limiting cases of this problem in which the equations reduce to simpler forms. Alternatively, you can learn to use Mathematica, Matlab or some other simulator to do it for you and explore the solution space.

Yeah, solving a coupled system of second order non-linear ODE is not covered in my textbook...

 

[bob012345]

Yeah, solving a coupled system of second order non-linear ODE is not covered in my textbook...

I think you can turn the second order equations into first order equations of different variables.

See if this numerical treatment makes any sense for you;

http://www.wiu.edu/cas/mathematics_and_philosophy/graduate/main.pdf

You will have to look up the unfamiliar terms such as the ERK method but that shouldn't be too difficult.