Numerically Simulating the Solar System in Mathematica

—The planetary motion within our solar system is a topic that has been studied for hundreds of years and has given rise to the science of astronomy. It is very important to know the positions of the planets in our solar system, as many of our current scientific research depends on it. Space exploration, for example, is a perfect example of when we need to know the exact positions of the planets in our solar system. Since it takes many years to send a rover or satellite to a planet, we will need to be able to predict the position of that planet many years into the future. Therefore, I present a second order Runge-Kutta simulation to predict the future position and velocity of the planets in our solar system based on Newtonian laws of motion. The equations of motion are implemented into a Mathematia script which animates the motion of each planet by generating a single static plot at each iteration within the while loop, stepping forward in time, re-plotting overtop the previous frame. This step-by-step numerical simulation is typically overlooked as an animation technique available in Mathematica. I herein provide an introduction to the software, an intuitive comparison of numerical vs analytical solutions to differential equations, and finally present the results of the simulation.


I. MATHEMATICA
Mathematica is a mathematical computation software that uses the Wolfram programming language, and it offers better symbolic manipulation than many other programming languages. Mathematica documents are called 'notebooks' and are organized into cells that can be individually evaluated. The Wolfram language has many built in functions, but it also allows the creation of custom functions. Functions make use of the fact that Wolfram is case-sensitive by using capital letters for built-in functions, which is why it is good practice to start custom functions and variables with lowercase letters. In order to define a custom function, variables is global across all notebooks. This works by assigning a serial number $ to the end of all variable names, making them unique.
Due to its user-friendly notebook format, Mathematica is commonly used as an elaborate graphing calculator and sometimes overlooked as a programming language for numerical simulations. Furthermore, Mathematica's selection of helpful built-in functions such as "Animate", "Manipulate", and "Dynamic", which allow the user to control input parameters (via a slide bar) for graphical plots of analytic functions, often cause users to overlook Mathematica's ability to animate numerical solutions in real time. The code presented in this paper shows a creative method to achieve step-by-step, real time animations by updating a particle's position in a while loop. To highlight the difference between analytical and numerical solutions, this paper begins with a discussion on differential equations, using the spring-mass system as an example. The spring-mass problem is solved both analytically and numerically to demonstrate the difference between the two approaches.

II. DIFFERENTIAL EQUATIONS: ANALYTIC VS NUMERICAL SOLUTIONS
A differential equation is any equation that involves the derivative of a function, with the derivative being the instantaneous rate of change of a function. Take for example: Notice the solution to the above differential equation is itself a function, ( ). This is true for all differential equations. Due to the time-derivative relationship between acceleration, velocity, and position, any physical system wherein the position is affected by the velocity or acceleration (i.e. essentially physical systems) will be described by a differential equation.
There are two approaches to solving differential equations, analytical and numerical. The difference between analytical and numerical solutions is that analytical solutions are more accurate, as their solutions are continuous functions, while numerical solutions are discrete functions. However, many real-life situations must be solved numerically, as their analytical solution either does not exist or is too difficult to solve for. A common real-world example where the analytical solution of a differential equation does not exist would the Blasius equation, used in fluid mechanics: ′′′ + 1 2 ′′ = 0 [1].
In order to demonstrate the two ways of solving differential equations, a simple spring-mass system will be used. Since a spring-mass system is much simpler than the solar system, it serves as an effective example to highlight the difference between the two methods.
The example system that will be used is a horizontal spring with a mass on the end of it, resting on a frictionless surface with no gravity involved. The goal is to solve for the spring's exact motion in time ( ), given an initial position ( = 0) = 0 and an initial velocity ( = 0) = 0 (schematic shown in Fig. 1). To solve for ( ) analytically, it is important to note that velocity is the derivative of position as a function of time, and acceleration is the derivative of velocity as a function of time. Therefore, = and = [2]. Thus, acceleration is equal to the double derivative of position as a function of time, = 2 2 . Newton's second law of motion states that = , and Hooke's law states that = − , where k is the force constant of a spring and x is the displacement of the spring from its equilibrium position. Setting the forces of the equations above equal to each other gives the equation − = . Setting acceleration to be the second derivative of position with respect to time, and dividing by , gives the following differential equation: Although it looks complicated, the equation simply means that the second derivative of (t) is the same function, here, Euler's identity has allowed to the exponential to be rewritten in terms of oscillating sine and cosine functions [3]. The unknown complex constants A and B are easily solved for when the initial conditions of the system are known. One can substitute ( ) → 0 and → 0, as well as ( ) → 0 and → 0 to yield two equations with two unknowns, solve for A and B, and obtain an exact solution for ( ) at all times. To highlight the difference between analytical and numerical solutions, the numerical solution to the same problem is now presented. Again, the goal is to find an unknown ( ) given a set of initial conditions: ( = 0) = 0 and ( = 0) = 0. In other words, "what is the motion of a mass after an initial perturbation?" First, the horizontal axis (time) is broken into discrete segments, as shown above in Fig. 2, and each ( ) must have some value at each point: { 1 , 2 , 3 …}. Next, it is assumed that the slope (i.e. derivative) between two neighboring points and +1 is simply the slope of a line that connects the two points, these are shown as red lines in Fig. 2. The slope of a line is given by mnemonic "rise over run", or slope ≈ +1 − +1 − . Since velocity is the derivative (i.e. slope) of position with respect to time, the following equation is obtained: Since acceleration is the derivative of velocity with respect to time: where Δ = +1 − has been substituted in. From Newton's second law, the acceleration of the mass must be -: Solving for +1 and +1 results in +1 = Δ + and +1 = − Δ + . These final two equations are well suited for computation. A while loop can be implemented to continuously calculate the new position and velocity { +1 , +1 } using the previous position and velocity { , }.
As simple as Euler's method is, it isn't the most accurate, as one side of the equation is centered with time, and the other side isn't. For example, in +1 − Δ = , the left side is centered halfway between +1 and , while the right side is centered at . Runge-Kutta is one solution to the issue, which attempts to solve the issue by centering both sides of the equation at the same through a half step using Euler's method to get and at +1/2 , which can later be used to more accurately solve for the next step. Applying this to the spring-mass system results in ℎ = − * Δ ℎ + and ℎ = * Δ ℎ + , which are then used to solve for the next step with +1 = ℎ * Δ + and +1 = − ℎ * Δ + .

III. NUMERICAL SIMULATION OF THE SOLAR SYSTEM
Moving on to the solar system, the mathematics is very similar to the mathematics in the spring-mass system. As there are many factors to consider in a system with multiple items, the system will be solved for numerically. Starting with Newton's second law of motion = , and Newton's universal law of gravitation = − 2 , the two equations can be substituted, then solved for acceleration, giving the equation: Noting that the mass of the planet cancels out. Using the fact that velocity is the derivative of position over time, the following formula is given: +1 − = , where can be replaced by to find the vertical vector. This equation can be solved for +1 , resulting in the equation: +1 = * + Now using the fact that acceleration is the derivative of velocity over time, gives the equation − 2 = +1 + , and solving for +1 gives the equation +1 = 2 * + . However, as the simulation is twodimensional, the equation for velocity needs to be split into and . To do this, the scalar is multiplied by its unit vector, ̂, to get vector ⃗. The unit vector is found by summing the and vectors of the planet, and dividing by the magnitude of ⃗ | | = √ 2 + 2 means that ̂= −̂−√ 2 + 2 . Together, this gives the equations: Using Runge-Kutta, the half steps of the previous equations are then taken, giving the final equations: where can be replaced by either or depending on the component that is being found. These last four equations are good for computational purposes, so given the initial positions and velocities of planets, as well as universal gravitational constant , the equations can predict the future positions the planets. It should be noted that there are factors of consideration that have not been accounted for; however, these factors do not produce a significant impact on the results. One assumption that is made about planetary movement is that the planets move in two dimensions (snapshots of the animation are shown in Fig. 3). In reality, the solar system exists in three dimensions, but due to the common orbital angular moment of the gas cloud that evolved into the solar system, the planets rotate in a single 2-D plane to a good approximation. At most, the orbital inclinations of the celestial planets only differ by 7 o with respect to on another, meaning very little accuracy is sacrificed in making this assumption [4]. A second assumption that is made is that the planets only feel a gravitational attraction to the sun. However, since the mass of the sun is far greater than any other object in our solar system, with the second most massive thing, Jupiter, being 1000 times less massive, a largest gravitational force a planet feels will be coming from the sun.
Before entering the initial values of all the planets into the code, the units of these values must be converted so Mathematica can run the simulation without having to deal with numbers that tens of digits long, thousands of times per second. To do this, all the distances were originally in meters, and were brought down by 10 9 . Times, originally in units of seconds, were divided by 33554.5, in order to set the velocity of the Earth at unity. The units of mass were kept at kilograms, as only the mass of the sun is ever used. The gravitational constant (units of m 3 kg -1 s -2 ) was converted to a value of = 7.514 * 10 −29 so as to incorporate the distance and time conversions used.

IV. CONCLUSION
Since the employed algorithm only used second-order Runge-Kutta and simplified the assumptions used (circular orbits in a two-dimensional plane), the presented simulation is less accurate than state of the art numerical simulations of the solar system, e.g. Runge-Kutta-Nystrom algorithms used by NASA [5]. However, the final script (attached at bottom of this paper), maintained impressive accuracy, matching true orbital patterns to within roughly ± 0.15 earth years for 40 earth years into the simulation. Fig. 3 shows the position of the planets from 0 to 100 years into the future. The overall accuracy of each planet is inversely proportional to its orbital period, meaning the planet with the shortest period, Mercury, is the first to deviate from its true orbit. If there exists even a slight discrepancy between the initial velocity and the initial position of the planet, then the planet will be thrown off axis very easily. This problem may be fixed if a higher order Runge-Kutta had been used. This simulation uses second-order Runge-Kutta, only taking the half step once, but a higher order, e.g. sixth-order Runge-Kutta could help even out the inaccuracies. Excluding Mercury, the other planets stay in orbit for ~500 years before being thrown off axis. The central purpose of this study, to develop a flexible method for simulating and animating physics in Mathematica was met. Although this paper only explicitly deals with the spring-mass system and celestial bodies orbiting the sun, the script (attached below) can be easily modified to simulate any number of physical phenomenon by substituting the necessary force in Eq. 6 and following the subsequent steps of Euler's (or Runge-Kutta's) method outlined in the paper.

ACKNOWLEDGMENT
The author would like to thank his family and friends, especially Kameron, for their help. He would also like to thank his high school physics teacher for piquing his interest in physics.