Computational rules might describe the evolution of the cosmos better than the dynamical equations of physics — but only if they are given a
“I think if he had had a more conventional education, if he’d come up through the ranks and had taken the standard physics courses and so on, maybe he would have done less interesting work.”
Today I’d like to share something I did as a project for a computational physics class. I had just learned how to solve ordinary differential equations (ODEs) in python, which gave me the tools to explore lots of physics questions: if I could describe the forces in a system using differential equations, I could set up my system in whatever configuration I want and see what happens when I let it go.
I decided to try this using the equations for Newtonian gravity between celestial bodies. I went the extra mile and learned how to animate my results, to better visualize how objects move through space. I also had the chance to try out some unusual configurations of objects, including ones that are chaotic and ones that use all three dimensions. Keep reading to see the equations and the neat animations they create!
One object
Let’s take a look at how gravity works on a single object first. We’ll write this as the acceleration due to gravity, or the second time derivative of our object’s position vector.
This is the simplest case, where we assume our object is orbiting around a stationary object at the origin with mass M. Our object’s acceleration always points toward the center, and gets stronger the closer to the center it is.
Now, in order to compute the motion of our object with that acceleration, we need to write it as a group of ordinary differential equations (ODEs). We’ll be using just two dimensions for now, so we’ll split up the equations for X and Y motion. We’ll need a summary of the object’s state at any given point: its x position, y position, x velocity, and y velocity. We need to provide the derivative of each of those state variables. Take a look and confirm that these all make sense given the acceleration above:
With this information, the ODE solver can work its magic: if we provide it with an initial state for our object to start at, it will integrate the above equations and show us how the state variables change. It’s basically a motion simulator!
In the following animation, our object is represented by the blue circle, and the black circle is the stationary mass at the center. I started my object away from the center, with a small upward velocity. The result is an elliptical orbit!
You could easily play around with the variables and see what changes. I set all the constants to unit values, but what if the mass at the center was greater? What if you start your object with a greater velocity? You might even see different types of trajectories under certain conditions. For example, if you go fast enough, the object will escape the system entirely!
But enough of that, let’s get real. Objects in real life don’t just stay stationary at the origin. They get accelerated by gravity too. Let’s see what happens when we let two objects move under each other’s gravity.
Two objects
The situation gets more complicated now that we need to keep track of the motion of two objects. We’ll call them objects 1 and 2, and their accelerations depend on the location of the other one:
It’s a lot to keep track of, but it’s not so complicated when we look at just the x and y directions. All we’re concerned with is the difference of the two objects’ coordinates. Here’s how it looks for object 1:
Using the same method as earlier, we generate a set of ODEs for every state variable, of which there are now 8 instead of 4. Instead of showing you all that, I’ll cut to the chase and show you the result when we set two objects of equal mass in motion.
When I give each object a small initial velocity in opposite directions, they dance around each other like so. They each have an elliptical orbit about their common center of gravity. You’ll see motion like this in the orbit of any two objects of similar mass, such as Pluto and Charon, or a binary star system. On the other hand, if we modify the variables so that one object is much more massive than the other, you’d end up with something that looks more like our first animation, but with the heavy object only moving slightly.
Three objects
Why stop at two? We can include as many objects as we want! We’ll just have to... re-write all our differential equations again. Yeah, you might be seeing now that this method isn’t sustainable for ever increasing numbers of objects. The acceleration of a single object in our system will now depend on the location of two others. It’s as simple as adding up the acceleration due to gravity from each object:
Or for a single component:
Furthermore, a third object means that motion can now happen outside the plane, so three dimensions are necessary to fully explore all the configurations. Three objects and three spatial coordinates means that my ODE solver needs a whopping 18 state variables, but it can handle the work just fine. I feed it all the differential equations and decide on some initial conditions to let it try.
How about we arrange three equal masses on the edge of a circle, equal distances apart? If we set them all off with the same radial speed, they should simply follow each other in a circular path, or maybe even dance around each other in elliptical paths like our previous animation. Here are the results of that:
Despite such a promising start, the motion ends up chaotic! It turns out that, for all but very specific initial conditions, three equal masses will never exhibit stable motion. My idea was sound: three objects positioned in a circle with the same radial speeds should result in simple orbit, but it’s an unstable equilibrium, and the slightest deviation (such as uncertainty in the initial conditions or the numerical solution) will result in chaos. If you want to read more about this phenomenon, the three-body problem is a well known topic in physics.
Don’t give up yet, there are still ways we can create a stable orbit using three objects. Think about the sun-earth-moon system: two small objects can orbit one large object with no problems. Choose one object to be our “planet”. Let’s increase the “sun” object’s mass by 100 times, and reduce the “moon” object’s mass by 100 times, then arrange them just so and let them go.
Another similar idea: have one small “planet” orbiting around two large binary “stars” that are orbiting each other.
With the right choice of masses, we can achieve stable orbits. We must only look to nature for examples.
Into the third dimension
We can get even crazier with our configurations of three objects if we leap out of the plane. Here’s one example called the Sitnikov Problem, which states that if two objects of equal mass orbit around each other, a third smaller mass can oscillate perpendicular to their plane of motion, threading their center of gravity like a needle.
Since we’ve already set up our ODE solver with three spatial dimensions, it’s not hard to choose the initial conditions that make this work. We’ll just need a little fancier plot to show it off.
If you think carefully about this configuration, you’ll probably notice that it’s an unstable equilibrium. The central mass must stay exactly on its line, or else it will be pulled away. Thankfully the computer can handle this just fine, since the central mass’s x and y coordinates can simply remain at zero.
The way I set it up, the two large masses are orbiting along a circular path. But if we change their speed so they orbit on elliptical paths, an interesting behavior arises: the central mass’s motion becomes chaotic. Still constrained to its line, but unpredictable. This is due to the way its gravitational potential changes as the large masses move in and out.
Pictured above is the position of the central mass over time. Indeed, testing confirms that its motion is sensitive to initial conditions, and thus fits the definition of chaos. Fun!
In conclusion...
My takeaway from this project is that differential equations are a very powerful tool for exploring the behavior of physical phenomena. Every simulation above is based on the differential equation for gravitational force, but you could create similar simulations for any dynamical system your heart desires.
Additionally, my increasingly complicated differential equations made me wonder: if I wanted to simulate a large number of objects interacting under gravity, what method would I use? What approximations would I make to simplify the calculations? Lots of interesting ideas there.
Today consisted mostly of writing flashcards for the math course I'm taking (exam next Friday! Eep) and programming. Currently working on an epidemiology-project, which is super interesting!
Didn't leave the apartment all day, so I'm really looking forward to getting out and meeting people tomorrow! 🌻
Current mood: wanting to be snarky towards trumpet players on my computational physics homework even though it’s been years since I last was in orchestra, but not being snarky towards them because I want the professor to take interest in me for research
