Modeling the Solar System, Part 2: Setting the Planets in Motion
December 2012
Our goal is to build a computational model of the solar system. In Part 1, we discovered where to obtain the necessary data and how to transform it for our needs. In the second part, we will set the planets in motion and test the accuracy of our model. All we really need to model our solar system is Newton’s law of gravitation.
\[ma_{x} = -\frac{GMmx}{r_{3}}\]If we initially assume we can use that equation to calculate the acceleration at each time step, then we can model the movement of the planets. At each step, the positions of each object are updated by adding the velocity and the velocities are updated by adding the accelerations. It is really that simple.
\[x = x + v_{x}; y = y + v_{y}; z = z + v_{z}\] \[v_{x} = v_{x} + a_{x}; v_{y} = v_{y} + a_{y}; v_{z} = v_{z} + a_{z}\]The equation for the law of gravity has a few important pieces: \(M\) is mass of the body not being updated at this step, \(m\) is the mass of the body being updated, \(r\) is the distance between the two bodies, and \(G\) is the gravitational constant.
We have not really discussed the mass yet. Again, to simplify the size of the numbers being used, we will use \(10^{24}\) kilograms as our unit. In these units, the mass of the Earth is 5.9736 and the mass of the Sun is 1,989,100 .
The only remaining issue is defining \(G\), the gravitational constant. The exact value of the constant is unknown, but it has been experimentally measured to a high degree of accuracy. We will use the tiny value
\[6.6738480 \times 10^{-11} m^{3} kg^{-1} s^{-2}\]Notice that the number is defined in terms of meters, kilograms, and seconds. Those are not the units being used in our model. We will need to redefine \(G\) in terms of AU, \(10^{24}\) kilograms, and minutes. This can be a little tricky if you are not careful—my original calculation sent the Earth hurtling through space faster than the speed of light. Even an error by one order of magnitude can cause the model to collapse.
I would encourage the interested reader to check the conversion, but the value used in the model is
\[7.17633288458 \times 10^{-17} AU^{3} 10^{24} kg^{-1} minutes^{-2}\]All units have been set and all variables have been defined. Just to be explicit, the update equation for acceleration is shown again. For a particular body, the result is obtained by summing over its interactions with every other body in the system. Only the update equation for the x component of acceleration is shown, but the other two components have identical formulations.
\[a_{x} = \sum_{i} - \frac{GM_{i} x}{ r_{3} }\]Testing the Model
To test the system, we load it up with initial data from February 9th, 1986, as we discussed in Part 1. Our model consists of the Sun, all 8 planets, the moon, and Halley’s comet. Through informal testing, this simple model of the solar system performs remarkably well. Even after several years the planets all seem to be in stable orbits and the moon continues to revolve around the Earth.
For a more definitive test, we have also collected data ten years after the initial set of data. How close is our model after simulating the solar system for 10 years—3,652 days according to this calculator? Not quite as close as I had hoped. The average distance of each body from its true location is 0.23 AU or approximately 30 million kilometers. That sounds like a huge error, but lets look a little closer.
Neptune is the largest offender, off by 0.59 AU, but that seems less egregious when you realize Neptune is 30 AU from the Sun. The planet that truly is the most difficult is Mercury. While its orbit only takes it about 0.35 AU from the Sun, our model is off by 0.21 after only 10 years.
The most surprising result is Halley’s comet. In the 10 years since February 9th, 1986, it goes from little more than half an AU from the Sun to over 20 AU. Our model estimates its position to within 0.06 AU.
The original goal was to have the model accurately predict one orbit of Halley’s comet. It is projected to again reach its nearest point with the Sun in 2061. At this point, the moon has drifted from the Earth and our prediction of the position of Halley’s comet is off by over 3 AU. Either the Earth will lose its only natural satellite within my life time or there is significant errors in the model. While my ego would like for the former to be true, I think mankind can rest assured it is the latter.
Sources of Error
Obviously the model is not perfect. There are likely many sources or error; here are just a few.
- Numerical rounding errors. Floating point arithmetic is not perfect and small errors can add up. We are modeling huge objects and vast distances, but a high degree of precision is still required. This can be difficult to handle.
- Initial velocity estimates. Our estimate of velocity is based on the difference between an objects place at two points in time, however, this is not the true instantaneous velocity. I imagine even a small error in initial velocity can lead to large errors given enough time.
- Timestep. A one minute timestep may be too coarse. The Earth moves nearly 2000 kilometers in that time. When objects are near each other, smaller timesteps may be necessary.
- Additional bodies. Our model currently only includes most of the major bodies in the solar system. There are thousands of other objects in the system, their absence from the model must introduce some amount of error.
Final Thoughts
I enjoyed translating the simple equations outlined in Feynman’s lecture notes into a working model. I have uploaded the data and C++ source code for anyone who is interested. Hopefully the code is understandable; I made an effort to comment the code and write it in a straightforward manner.
Further extensions could add a graphical component to visualize the orbits or focus on improving the accuracy of the model. I think something like this would make a great project for a programming inclined science student. If anyone has worked on something similar or has ideas for improving the model presented here, I would love to hear about it.