We will solve the equation of orbital motion, which is a second order ordinary differential equation, using the second order Runge-Kutta method. We may write this equation as a vector equation:
where bold faced r is a vector:
and non-boldfaced r is the magnitude of the vector (i.e.
the distance between planet and Sun in this case). Note that this
means that r/r is a unit vector pointing
from Sun to planet.
Note that one of the strengths of MATLAB is that it can handle
vector operations very easily. This can be used to achieve some
simplification in the writing of the program if you like.
For simplicity, we will carry out the calculation using
length units of Astronomical Units and time units of Years. This has the
virtue of simple constants and solutions that are written using small
numbers. In this system of units,
To solve this equation, we need to supply initial conditions. Let the calculation start at x=1 AU, y=0 AU, and z=0 AU. Let the initial velocity in the x and z directions be 0, and write a script that solves this differential equation for an initial velocity in the y direction that is input by the user.
Please make a subdirectory for this exercise under your
directory ~schloerb/ph281/username and then copy the
following into that subdirectory:
Calculate a family of 3 orbits, starting at the same initial position but
with different velocities. One of the velocities should be that which gives
a circular orbit. The other two velocities should be selected to be a bit
below and a bit above the circular velocity. Make a graph with all three
orbits. How does the perihelion position change?
Extra Credit