Mathematica: Differential equations for double pendulum

[furrer]

Hello

Im new here, I hope I'm not disturbing anyone.

Following this guide below, I am trying to find two 2. order differential equations, one for q1'' and one for q2'', describing the movement of the double pendulum. I have no problems with the mathematics, but when the guide tells me to use Mathematica, I have some problems in getting the same result as in the guide and elsewhere.
Could you help me?
I'm attatching my recent try in mathematica.
Note: I'm not trying to numerically solve the differential equations, as done later in the mentioned guide, I'm just trying to "create" those differential equations.

The guide mentioned is:
http://myphysicslab.com/dbl_pendulum.html

Sorry for my bad english.

 

[Simon_Tyler]

So, as the guide says, you need their equations 1-4 and 13,16

eq1 = x1''[t] == -q1'[t]^2 L1 Sin[q1[t]] + q1''[t] L1 Cos[q1[t]];
eq2 = y1''[t] == q1'[t]^2 L1 Cos[q1[t]] + q1''[t] L1 Sin[q1[t]];
eq3 = x2''[t] ==  x1''[t] - q2'[t]^2 L2 Sin[q2[t]] + q2''[t] L2 Cos[q2[t]];
eq4 = y2''[t] == y1''[t] + q2'[t]^2 L2 Cos[q2[t]] + q2''[t] L2 Sin[q2[t]];
eq13 = Sin[q1[t]] (m1 y1''[t] + m2 y2''[t] + m2 g + m1 g) == -Cos[q1[t]] (m1 x1''[t] + m2 x2''[t]);
eq16 = Sin[q2[t]] (m2 y2''[t] + m2 g) == -Cos[q2[t]] (m2 x2''[t]);

Then, I think the easiest way to solve the system is to use substitute equations 1 to 4 into equations 13 and 16, then solve the result.

eq13 && eq16 //. ToRules[eq1 && eq2 && eq3 && eq4] // Factor;
soln = Solve[%, {q1''[t], q2''[t]}];

The solutions require a bit of cleaning. I found it easiest to reproduce the results in the guide by dealing with the numerator and denominator separately. The following code reproduces their results exactly.

q1''[t] /. soln[[1]] // Together[#, Trig -> True] &;
q1num = Collect[FullSimplify[Numerator[%]], {Sin[_]}, Factor]
q1den =Factor[Denominator[%%]] 

q2''[t] /. soln[[1]] // Together[#, Trig -> True] &;
q2num = Collect[FullSimplify[Numerator[%]], {Sin[_]}, Factor]
q2den = Factor[Denominator[%%]]

I wrote a very quick and dirty bit of helper code to turn the variables used above that are nice in forums into ones that look nice (e.g., m1 -> m₁)

symRep[expr_] := expr //. s_Symbol?(StringMatchQ[SymbolName[#], 
       WordCharacter ~~ DigitCharacter] &) :> (
      Subscript[Symbol[#1 /. "q" -> "\[Theta]"], #2] & @@ Characters[SymbolName[s]])

Using it and TraditionalForm to display the above results gives

q1''[t] == q1num/q1den // symRep // TraditionalForm
q2''[t] == q2num/q2den // symRep // TraditionalForm

Note, you also might want to look at
http://demonstrations.wolfram.com/DoublePendulum/
as well as the derivation of the differential equations using the Euler-Lagrange equations
http://scienceworld.wolfram.com/physics/DoublePendulum.html
http://en.wikipedia.org/wiki/Double_pendulum#Lagrangian

 

[furrer]

Thank you Simon, I'ill look at it a bit later, right now I'm not at a pc with Mathematica. It looks perfect though, thank you!

As you probably could see by my coding, I'm very new to Mathematica.
I knew about the Euler-Lagrange method, but that's above my level in physics.
I'm using this for a final high-school paper in physics and mathematics. I know the double pendulum is a very (and too) complicated system for a high-school-student, but when I found that Newtonian method, I thought it might just go right, and I really did only pick my subject (non-linear equations and chaos) for the mathematics!

Maybe I'll be back with more questions later, if it's ok.

 

[Simon_Tyler]

Anything beyond very basic differential equations is beyond the average high school course. So you get my admiration for attacking this problem!

I didn't look at your code before, but I have now.
You actually got quite close! Your main error was using sin and cos instead of Sin and Cos.
In standard input and output, all built-in functions start with a capital letter.
(You might want to work through the "getting started" tutorial http://reference.wolfram.com/mathematica/tutorial/IntroductionOverview.html)

So your mistake prevented Mma from using any trig identities in the simplification.
I've attached a version of your notebook with the minimal correction.

 

[furrer]

Thank you. Soon I'll have two free weeks, where I am to write a paper about non-linear dynamics and chaos theory. With the level of education I have, I'm taking a illustrative route (with some comparing to fractals), but I thought I at least should find the equations of motions for the double pendulum, after which I could examine it's periodic and chaotic motions by simulations. Luckily, a lot of such simulations have already been programmed for Mathematica, so I'm avoiding having to write a simulation-program!

I didn't look at your code before, but I have now.
You actually got quite close! Your main error was using sin and cos instead of Sin and Cos.
In standard input and output, all built-in functions start with a capital letter.

Ahh of course. I had actually looked a lot at that tutorial, but I forgot about the capital letter. In school we are working with something as old as TI-InterActive (nearly older than me! - hehe), where you don't have to do this. Silly me.

Thank you for the attachment and your earlier solution, which was great, also in understanding how to do the same for other problems. This help learned me a lot about using Mathematica, so thank you for this too.

Do you know, that when typing NDSolve, what method and what stepsize Mathematica is using? I know it's possible to command what method and stepzise to use, but I'm a bit curious, and I couldn't find it anywhere in the Virtual Book.

Lastly...
In my paper I would like to examine the fractal nature of the double pendulum movement, again avoiding energy, as this requires Euler-Lagrange. I have an idea about a graph, with θ1(start) as x and θ2(start) as y. Than for specific L1, L2, g, m1, m2 I would like mathematica to plot a white pixel if θ1 is positive after for example 20 s and a black pixel if θ1 is negative (using 4th order Runge-Kutta). I'm afraid this is far beyond my abillities when it comes to programming, but maybe moduls like this have already been made?

Thank you for the help, I hope it's ok that I'm asking more questions.

 

[Simon_Tyler]

Hi Furrer,

Sorry for the following link dump - I'm a bit short of time today.

NDSolve uses various heuristics to choose an appropriate method of numerical integration and is able to dynamically change method if difficult parts are encountered in the integration. It's loosely described at http://reference.wolfram.com/mathematica/tutorial/NDSolveIntroductoryTutorial.html

For your last point, what you want is probably something like a Poincare Section.
The way this is normally done is to use the StepMonitor option or EventLocator method of NDSolve to create the data for the section.
https://groups.google.com/forum/#!topic/comp.soft-sys.math.mathematica/ZWGPZMnW5Fo
and there's a good example in
http://reference.wolfram.com/mathematica/tutorial/NDSolveEventLocator.html
You can also see some working code at
http://demonstrations.wolfram.com/search.html?query=Poincare section
There is also EquationTrekker that can do Poincare sections (although, I haven't used this package before) http://reference.wolfram.com/mathematica/EquationTrekker/tutorial/EquationTrekker.html

Another way to examine the chaos in the double pendulum is to see where m2 ends up after a given time and small variations of the initial conditions. I think I saw a post about this on http://blog.matthen.com, but the servers are apparently overloaded at the moment! You can view the pages via google cache if you want. He placed some code here: http://pastebin.com/ndFf0yXm

 

[furrer]

Hi Simon

No problem, links are great!

Ok, thank you for the explanation and the link, which will be further investigated.

Actually I wasn't thinking about a Poincare Section, which is far above my level. As far as I know, this requires finding the energy in the system, but I might be wrong. As I avoided the Euler-Lagrange, I think this might not be possible? But I will look at the links later and reply more concretely (the time difference between me (Denmark) and you is just brilliant!). After all, it would be great, if it was possible, as the Poincare Section is very interesting and in my opinion near the essence of chaos!

I was thinking about this instead:
https://freddie.witherden.org/tools/doublependulum/report.pdf (page 5)
It might not be a real fractal, but I think it could be compared to a fractal? It just seems hard to code such an application in Mathematica, and he didn't show all his coding in this paper.
And this:
http://tabitha.phas.ubc.ca/wiki/images/archive/3/37/20080811183757!Double.pdf
Which my teachers luckily say, I'm allowed to investigate without coding myself (i.e., just use their "pictures"). (as this one is far to complicated to code; and their code is not to be found I think).
But I would really prefare to be able to also make some simulations like the one from the "freddie.witherden"-homepage, as this gives me more freedom and a more independent paper.

Your last idea have I already considered, the butterfly effect is always interesting to investigate. The little modul was good, but I also found this one myself:
http://demonstrations.wolfram.com/DoublePendulum/
Both might be good to use. And you are talking about q2, right?

Again, thank you for the good response, I'll reply again later, considering the links about the Poincare Section of the double pendulum.

2nd reply:
"Because this is a Hamilton non-dissipative system, inital energy of the system is a constant of time and initial energy is a function of initial conditions. To get a real poincare diagram i must repeat the procedure described above for different initial conditions, but for the same energy level."
This quote is from here: https://groups.google.com/forum/#!topic/comp.soft-sys.math.mathematica/ZWGPZMnW5Fo
After all, this is what I'm trying to avoid!
So I would really hope, that you have some ideas how to implement my other idea in Mathematica!

Thank you for help and understanding kindness.