Online study group for dynamical systems theory and innovative trajectories

[mmeijeri]

Anyone interested in an online study group for this?

I am.

This thread is for a study group for the following online book: Dynamical Systems, the Three-Body Problem, and Space Mission Design by Koon, Lo, Marsden, and Ross. Discussion of the mathematics and physics involved as well as the applications to alternative exploration architectures are on topic for this thread.

[Hop_David]

Anyone interested in an online study group for this?

I am.

This thread is for a study group for the following online book: Dynamical Systems, the Three-Body Problem, and Space Mission Design by Koon, Lo, Marsden, and Ross. Discussion of the mathematics and physics involved as well as the applications to alternative exploration architectures are on topic for this thread.

I'm printing out pages 1 through 56 (what Adobe Acrobat calls pages 15 through 70 of the pdf). This is the first three chapters. I will read until I reach regions opaque to me and then ask questions.

[engstudent]

Interesting.

[Hop_David]

Hit my first speed bump on page 8. Don't know how to replicate KoLoMaRo's notation in this venue so I'll define my own notation:

x' = dx/dt
y' = dy/dt
p = position vector = (x, y)
p' = v = dp/dt

Now I had thought kinetic energy was 1/2 * m * v^2

where v = |v| = sqrt(x'² + y'²)

This would give kinetic energy as
1/2 * m * (sqrt(x'² + y'²))²
which is
1/2 * m * x'² + y'²

However they give kinetic energy as
1/2[(x' - y)² + (y' + x)²]

The - y term and the + x term mystify me.

Even though the equations are somewhat mysterious I think I understand "Thus at a given energy level E, the tyhird body can only move in the region given by the inequality E - U >or= 0; this is called the Hill's region and is obtained by intersecting the graph of the effective potential with a horizontal plane."

Figure 1.2.3 illustrating this sentence is really neat.

I am reading ahead and studying illustrations to try absorb what I can. I am hoping the opaque equations will become clear as I study them. I am pleased to add lyapunov orbits, lissajous orbits and Conley-McGehee tubes to my vocabulary.

[mmeijeri]

Now I had thought kinetic energy was 1/2 * m * v^2

They are using a rotating frame with origin at the center of mass. They are also assuming the primaries move in perfectly circular orbits around the center of mass, as opposed to the more general case of elliptical orbits. Let A(t) be the 2x2 matrix defining the transformation from the moving frame with respect to an inertial reference system in the joint orbital plane with the origin at the center of mass, then we have:

A(t) = R(omega * t),

with R(theta) the 2x2 rotation matrix for a counterclockwise rotation through an angle theta and omega the angular velocity of the rotating frame. The rotating coordinates can then be described as a function of the inertial coordinates as follows:

r^~ = A * r

with r being the position vector in the inertial system and r^~ the one in the rotating frame.

Equivalently,

r = A^-1 * r^~

Differentiating we get

r = A^-1' * r^~ + A^-1 * r^~'

by the product rule for differentiation. The first term on the left gives you the additional terms.

In coordinates:

A-1 =

cos omega * t -sin omega * t -sin omega * t -cos omega * t

A-1' =

omega * -sin omega * t -cos omega * t -cos omega * t +sin omega * t

Geometrically, the first column vector of A describes the rotation of the x-axis. Differentiating, we get the tangent at a point on a circle, which is precisely perpendicular to position vector of that point. Substituting into the formula for r~, we get

r' = omega * r^~_x * R(-omega * t) * (0 1)^T - omega * r^~_y * R(-omega * t) * (1 0)^T + R(-omega * t) * v^~

or

r' = omega * R(-omega * t) * (-r^~_y +r^~_x)^T  + R(-omega * t) * v^~ = R(-omega * t) * (omega * (-r^~_y +r^~_x)^T + v^~)

This gives the following expression for the kinetic energy:

E_kin = 1/2 * m * (r',r') = 1/2 * m * (R(-omega * t) * (omega * (-r^~_y +r^~_x)^T + v^~),R(-omega * t) * (omega * (-r^~_y +r^~_x)^T + v^~))

Since rotation preserve angles, we can leave out the rotation matrices:

E_kin = 1/2 * m * (omega * (-r^~_y +r^~_x)^T + v^~),omega * (-r^~_y +r^~_x)^T + v^~) = 1/2 * m * omega² ((-r^~_y +r^~_x)^T + v^~,(-r^~_y +r^~_x)^T + v^~)

Units of time have been chosen in such a way that omega = 1, which leaves us with the formula given by KoLoMaRo.

[mmeijeri]

Yikes, typing those formulas into the editor is a lot of work!

[Hop_David]

Yikes, typing those formulas into the editor is a lot of work!

Thank you! I hope you will be patient with me while I study your post for awhile.

[mmeijeri]

Oh wait, there's a much easier way of seeing it: the extra term isn't just caused by the rotation of the rotating frame with unit angular velocity, it is the velocity of the rotation. It's perpendicular to the position vector, points in the right direction and its length is proportional to the length of the position vector.

[Hop_David]

I've used rotation matrices. And I know the product rule for derivatives. But I don't recall differentiating a rotation matrix. My first thought was angular velocity is constant so A' would be zero.

But I can see differentiating each trigonometric expression within a 2x2 rotation matrix would give a new rotation matrix 90 degrees from the original. And that, with constant angular velocity and circular orbits, the velocity vector would be proportional and perpendicular to the position vector.

Your next reply pulled away the veil. I can see in a rotating frame that there would be an apparent motion perpendicular and proportional to the position vector from barycenter.

I expect angular velocity will come up often. How about calling it "w" rather than "omega"? that would save some keystrokes.

[Downix]

I wish I had a head for this.  Instead I am good with boolean and electronics systems.  But it is a very good read none the less.

[Hop_David]

I put up a CR3BP model
(based on Bob Jenkin's Java Module orb.jar orb.jar source)

Or this might be called CR2+10BP since I have two massive bodies and 10 massless objects.

The mass of the central body is 1000 and the less massive body is 250.

So
µ = .2
1 - µ = .8

(does µ show as mu? I am going to try to use some non-ASCII characters. Holler if they come out as gobbledygook)

At the bottom of page 9 & top of 10: "Thus, at a given energy level E, the third body can only move in the region given by the inequality E - UŻ ≥ 0"

Here's a quick illustration showing parts of the Hill region that I'd imagine to have close to the same effective potential as L1:


The massless orbits in my CR3BP simulation don't remain in the Hill Region. They either drop to olive shaped orbits about the sun or they orbit about Jupiter making a big doughnut. They seem to drift into either the Sun or Jupiter realm. (as the authors call regions in Figure 1.2.3)

I know an ordinary elliptical orbit conserves energy as potential energy waxes while kinetic wanes (and vice versa).

The orbits that fall towards the sun don't seem to have constant angular momentum. For a time both the apohelions and perihelions seem to shrink as Jupiter comes from behind during every 5th (or so) apohelion and steals velocity. Then after a time the massless bodies come up from behind Jupiter as they ascend towards apohelion and then Jupiter seems to boost their angular momentum.

I would guess the angular momentum  of the entire system remains constant as the tiny bodies borrow or lend neglible amounts of angular momentum from the sun or Jupiter.

But the sentence bottom page 9, top page 10 stills bothers me. I wouldn't think the E of a body changes much as it drifts from L1 into the sun's realm. But to my eyes it sure looks like it leaves the Hill Region.

[Hop_David]

I wish I had a head for this.  Instead I am good with boolean and electronics systems.  But it is a very good read none the less.

How do you think I feel? I graduated in the upper 60% of my high school class and majored in art during my short time at a university. An obsession with M.C. Escher led me to a math hobby. I've learned much of the math and physics I know from the internet, library books and books purchased at yard sales. I often feel like a non-swimmer in the middle of a deep lake.

Rummaging in my shed for helpful books I discovered Vector Calculus third edition by Jerrold E. Marsden and Anthony J. Tromba. Jerrold Marsden is one of the co-authors of the pdf Mmeijeri has linked to! I wonder if Jerrold is related to Brian Marsden of the Minor Planet Center.

[mmeijeri]

I put up a CR3BP model

I'm getting an error loading that page, I'll see if it has something to do with my Java settings.

The massless orbits in my CR3BP simulation don't remain in the Hill Region. They either drop to olive shaped orbits about the sun or they orbit about Jupiter making a big doughnut. They seem to drift into either the Sun or Jupiter realm. (as the authors call regions in Figure 1.2.3)

I don't think angular momentum of the infinitesimal particle can be constant since it is not moving in a central force field. Total energy and(?) angular momentum of the full three body system should be conserved.

Bear in mind that your numerical approximation may not conserve either angular momentum or energy. Explicit Euler for instance increases both energy and angular momentum, and I believe implicit Euler does the opposite. An average of the two might work well, but there is a better way called Euler-Cromer. The trick is to update the positions first and then use the updated positions in the calculation of the forces and use that for updating the velocities:

x_n+1 = x_n + v_x * dt
y_n+1 = y_n + v_y * dt
a_x_n+1 = F_x(x_n+1,y_n+1)/m
a_y_n+1 = F_y(x_n+1,y_n+1)/m
v_x_n+1 = v_x_n + a_x_n+1 * dt
v_y_n+1 = v_y_n + a_y_n+1 * dt

If you look closely, it is an implicit method, but because you do not need to know v_x_n+1 and v_y_n+1 for the position updates you can still solve it explicitly. That's why it's called semi-explicit.

Euler Cromer exactly conserves angular momentum. The energy is not exactly conserved, but oscillates with the same period as the orbit. If you take your timesteps small enough, the amplitude of the oscillation can be made arbitrarily small. For small numbers of orbits RK4 should also work well, but this does not conserve angular momentum in the long run. This matters if you want to do long-term simulations of the solar system.

Euler Cromer is a simple example of a so-called symplectic integrator, and there is a whole body of mathematics behind it, of which I know only very little. Unsurprisingly, if you google for this you'll find articles by Marsden, who is a giant in the field of differential equations and dynamical systems. I think we should stick with using Euler-Cromer for now until we run into trouble. It's nice and simple.

I'n the spreadsheet I'm working on now I calculated the orbit of the moon this way. I'm deliberately forcing myself to stick to Excel and VBA because I want to get more familiar with them. Normally I would do this sort of thing in C++. To keep things easy I started with explicit Euler, and as I had expected it conserved neither energy nor angular momentum. Euler Cromer worked spectacularly well, especially after I had fixed a few incorrect minus signs, which had led to the slightly implausible prediction of 21.6 day for the period of the moon. I think at some point I'll be forced to switch to Java, because VBA is very limited.

[mmeijeri]

David, could you put the required classes in your web directory? I think they may be in your class path on your own system, but I can't find them from here.

[Downix]

I wish I had a head for this.  Instead I am good with boolean and electronics systems.  But it is a very good read none the less.

How do you think I feel? I graduated in the upper 60% of my high school class and majored in art during my short time at a university. An obsession with M.C. Escher led me to a math hobby. I've learned much of the math and physics I know from the internet, library books and books purchased at yard sales. I often feel like a non-swimmer in the middle of a deep lake.

Rummaging in my shed for helpful books I discovered Vector Calculus third edition by Jerrold E. Marsden and Anthony J. Tromba. Jerrold Marsden is one of the co-authors of the pdf Mmeijeri has linked to! I wonder if Jerrold is related to Brian Marsden of the Minor Planet Center.

Your story sounds similar, save I was at the bottom of my class, not due to being dumb, but that my class was so stacked. (the top 20% all had 99% or higher)  Went to univ for a year for filmmaking, then got into video editing discussions with an Amiga fanatic.  Showed him some tricks, then got to discussing the hardware inside of it, and got me into studying electronic chips and how they worked.  Self tought, altho I am planning on going back for an EE degree now.

[Hop_David]

David, could you put the required classes in your web directory? I think they may be in your class path on your own system, but I can't find them from here.

Checked with my ftp program.   orb.jar is in the railroad folder on clowder.net/hop/railroad. Ran the page from my co-worker's computer and it worked OK. This is troubling. I was hoping to make sim pages accessible to any person browsing the internet. Well, here's a couple screenshots:

[simonbp]

You gotta let it evolve for a while... I got two (red and blue) in semi-stable orbits around the secondary body, and the rest ejected...

Simon

[Hop_David]

The trick is to update the positions first and then use the updated positions in the calculation of the forces and use that for updating the velocities:

x_n+1 = x_n + v_x * dt
y_n+1 = y_n + v_y * dt
a_x_n+1 = F_x(x_n+1,y_n+1)/m
a_y_n+1 = F_y(x_n+1,y_n+1)/m
v_x_n+1 = v_x_n + a_x_n+1 * dt
v_y_n+1 = v_y_n + a_y_n+1 * dt

I recall making a 2 body model with a spread sheet using the above. Instead of an ellipse, I got a slowly growing spiral. Which is sort of what I expected since I knew I was piecing together a multitude of small parabola fragments. Jorge Frank told me I was using a first order Euler integrator . He suggested I get the acceleration at the midpoint of the line segment connecting  (xn, yn) and (xn+1, yn+1). He said this was a second order Runge Kutta integrator.

Following Jorge's suggestion I was able to get somewhat better results. But I didn't get past huge, unwieldy spreadsheets that I had little confidence in.

Then a friend pointed me to Bob Jenkin's  web site. Bob seemed to have given orbit modeling some study. And using his orb.jar was a whole lot more convenient than my horrible spread sheets. Bob talks about multi step methods here.

I am still disappointed and scratching my head that the CR3BP model based on Bob's orb.jar didn't run for you.

This page has links to a variety of pdfs. In his Introduction to the N-Body Problem Jay James seems to be saying the system's angular momentum stays constant as well as it's energy. If you choose a convenient coordinate system, the system's center of mass stays put.

To be honest, I don't know how well Jenkins' model conserves energy or angular momentum. In his notes I saw no mention of Euler Cromer.

Googling VBA . . . Visual Basics for Applications? Terra incognita for me. I am comfortable with Excel, though.

Hopefully we will find some modeling methods that are accessible to everyone.

[simonbp]

He suggested I get the acceleration at the midpoint of the line segment connecting  (xn, yn) and (xn+1, yn+1). He said this was a second order Runge Kutta integrator.

What he suggested was a Symplectic integrator, specifically a Velocity Verlet. Symplectics are far better for orbital integration, as they conserve the total energy, unlike a Runge-Kutta. On the other hand, Runge-Kuttas (fourth order is usually a good compromise) are better when the total energy is changing, i.e. when a rocket is burning....

Simon

[mmeijeri]

I recall making a 2 body model with a spread sheet using the above. Instead of an ellipse, I got a slowly growing spiral. Which is sort of what I expected since I knew I was piecing together a multitude of small parabola fragments. Jorge Frank told me I was using a first order Euler integrator.

Yes, that's what you get with explicit Euler. But note that the scheme above is not explicit Euler, because you use the already updated positions in your velocity update. That's what makes the scheme symplectic. It's this change that makes it Euler-Cromer vs ordinary explicit Euler.