Simulation of the Solar System and Evaluation of Runge Kutta Methods

By admin On March 25th, 2010

Introduction

For this project, I created a simulation of the n-Body Problem and used this to simulate the Solar System. I extended this project to an analysis of second- and fourth-order Runge Kutta ODE solvers by verifying realistic outcomes and conservation laws.

Sun, Eight Planets, and Earth’s Moon (2nd Order)

I simulated the Sun, the eight planets, and Earth’s Moon using a second-order Runge Kutta ODE solver. I initialized the system by assuming that all the planets are in uniform circular motion about the center of the Sun (In reality, they are in almost circular motion and they do not quite orbit the center of the Sun, but this is a good approximation since the mass of the Sun is much greater than the masses of the planets). Using the two equations for uniform circular motion:

a = v2/r
v = 2 π r / T

where a is the acceleration due to the Sun’s gravity, v is the initial velocity, r is the radius of the planet’s orbit, and T is the period of the planet’s orbit. Using these two equations, I was able to solve for v in terms of a and r only. Since these are known, we was able to easily calculate the initial velocities of all the planets.

As for the moon, I similarly calculated the initial velocity of the moon about the Earth and (vectorially) added the Earth’s velocity to this in order to ensure that the moon continued to orbit the Earth as the Earth orbited the Sun.

At each time step, the position and velocity was computed by using Newton’s first law (F=ma) and Newton’s law of gravitation to produce:

x’ = v

v’ = G m / r2

A plot of the system is shown below (The time step used was 20 Earth days).

The moon indeed orbited the Earth, as shown below:

In addition, I verified that energy was conserved by calculating kinetic and potential energies at each iteration. Below, the green plots the total kinetic energy and blue plots the total potential energy. Red plots the sum of the kinetic and potential energies.

I also verified that the total linear momentum was conserved:

Finally, I verified that the total angular momentum was conserved:

Sun, Eight Planets, and Six Moons (2nd Order)

Next, I added one moon to each planet (except for Mars and Venus, which do not have moons). The initial velocities of all bodies were calculated in the same manner as described above.

A plot of the system is shown below:

As can be seen in the above image, none of the moons have stable orbits except for Luna (Earth’s moon) and Titan (Saturn’s moon). I looked to see if energy was conserved in this system. Below, the green plots the total kinetic energy and blue plots the total potential energy. Red plots the sum of the kinetic and potential energies.

Surprisingly, the total energy is conserved. This indicated that simply using a smaller time-step may solve this problem. However, using a time step as small as 2 hours improved the stability of the moons’ orbits only slightly.

I also verified that the total linear momentum was conserved:

Finally, I verified that the total angular momentum was conserved:

Sun, Eight Planets, and Six Moons (4th Order)

Upon further examination the Solar System, I determined that Luna and Titan have the largest orbits of all the moons. Hence, I reasoned that using a higher-order ODE solver might resolve the problem with instabilities in the orbits of the moons. For this reason, I implemented fourth-order Runge Kutta.

A plot of the system is shown below:

As can be seen in the above image, all of the moons have stable orbits except for Cordelia (Uranus’ moon). I noticed that the period of Cordelia’s orbit is an order of magnitude smaller than any other moon, indicating that either using a smaller time step or a higher-order ODE solver might resolve this problem. The time step I used to generate the above image was 1 Earth day — I tried using a time step of 2 hours, but the simulation did not complete after running it for three hours.

I looked to see if energy was conserved in this system. Below, the green plots the total kinetic energy and blue plots the total potential energy. Red plots the sum of the kinetic and potential energies.

Energy is conserved in this system. Since using a smaller time step may not be feasible, it may be necessary to use a more advanced ODE solver such as ode45.

I also verified that the total linear momentum was conserved:

Finally, I verified that the total angular momentum was conserved:

Fun Fact

If I plot the motion of the Sun in my simulation, it looks like this:

Note that the Sun appears to be orbiting about some point — and the period of this orbit is nearly identical to Jupiter’s (this is because Jupiter is to massive). When astronomers seek to identify which stars have planets orbiting them, they look to see which stars are wobbling. However, this method only works to identify Jupiter-sized planets. If intelligent life existed somewhere else in the universe and they used similar techniques to identify planets orbiting other stars, they wouldn’t be able to notice the Earth, since our planet is too small to cause a significant wobble in our Sun.

Conclusions

In conclusion, I found that both second- and fourth-order Runge Kutta ODE solvers conserve energy in the n-Body problems I considered — however, I noted that second-order Runge Kutta did not yield realistic results for most of the moons since the periods of their orbits are small. Fourth-order Runge Kutta yielded more realistic results, although it failed for Cordelia, which has a very small orbital period. The best way to improve my system would probably be to use a more sofisticated ODE solver such as ode45.

Code

2