Joel Berrier

Hope College

Research Professor - Dr. Peter Gonthier

Supported by the NSF-REU and NSF

For my research project this summer I continued the development of the computer program originally created by Michelle Ouellette during past summer research programs. The purpose of this computer program is to simulate the birth and distribution of pulsars in the galaxy starting with several randomly determined parameters in a cylindrical coordinate system. These parameters are randomly generated from distributions using information taken from the currently available pulsar literature. After these initial parameters are determined the program continues by evolving these stars through to their present state using a numerical integration method, specifically a Runge-Kutta calculation.

Over the course of the summer I hoped to begin the development of a useful 5^{th} order Runge-Kutta-Fehlberg method to replace the previously used 4^{th} order Runge-Kutta routine. The program uses a Runge-Kutta method to integrate the equations of motion while it is evolving the star to its current state from its initial, randomly assigned, conditions. This integration can cover a span of up to 10^{9} years. The 4^{th} order code also had a crude adaptive step size to help adjust itself to the current error in the energy.

This 5^{th} order Runge-Kutta method is known as an embedded Runge-Kutta. Fehlberg originally developed this method. Included in this embedded Runge-Kutta routine is a method to adapt the step size by using the truncation error between a 4^{th} order and 5^{th} order calculation of the relevant quantities. The relevant quantity in our program is the error in the energy. We used the difference between these 4^{th} and 5^{th} order errors in energy to determine our next step size.

We decided to attempt to use this method to create a more efficient program. We believed that this would be accomplished by taking a larger possible step size, thereby reducing the total integration time, and controlling the possible error in the energy. Also we decided to use the Cash-Karp parameters for the Runge-Kutta, due to their higher level of accuracy, instead of the parameters originally created and used by Fehlberg. Using this method we hoped to decrease the total integration time and maintain an acceptable level of error in the energy.

Thus far the quality of our results with the embedded 5^{th} order
Runge-Kutta-Fehlberg routine have been less then what we had desired. The 5^{th} order
Runge-Kutta code produces a far larger number of bad events then in the previously used
4^{th} order code. A bad event is defined, within the code, as an event in which the
fractional error in the energy has grown to a value greater then a set parameter, we currently
use a value of 1.0*10^{-5}. The previously used 4^{th} order routine also
currently runs approximately 2.5 times faster than the recent version of the 5^{th}
order code.

The adaptive step size is controlled by the difference between the fractional error
in the energy as calculated with both the 5^{th} and 4^{th} order methods. This
quantity is scaled to our problem by dividing it with the current energy from the 5^{th}
order calculation then dividing it again by our error controlling variable EPS. EPS is the
desired accuracy of our routine. The final product of these calculations, designated as
errmax in the code, is used in standard step size calculations given below. If the errmax
is greater then 1.0 the program begins to reduce the step size, but it reduces it to no less
then one-tenth of the original step size. If the value is greater than an error control value,
approximately 1.889*10^{-4} in our program, see calculation below, we increase
the step size by a small margin. If the value of errmax is less then the error control
value we raise the step size to five times the previous step size. The standard equations
for step size growth and step size reduction in the 5^{th} order embedded
Runge-Kutta-Fehlberg method is as follows: Growth ^{(-0.20)}, reduction
^{(-0.25)}.

The error control parameter is determined by the following formula: . The variable safety is used to prevent to rapid a change in the step size during both growth and reduction. This value needs to be just a few percent less then unity. For our problem we have used 0.9 as our safety value. The error control is used to prevent too much growth in step size as the errmax approaches 1.0.

The adaptive step size is not as responsive to the changes in error as we would like
it to be. We had hoped that it would adjust itself to a suitable size to control the growth of the
error. As the error in the energy increases the step size should shorten to prevent any large
growth of the error in the energy. The value of the step size during the integration, however,
selects a value that it settles on within the first few steps of the total integration. Then
the step size remains at that value for the rest of the integration. This portion of the code
is the major reason why we generate so many bad results, as well as why we take longer then
the fourth order Runge-Kutta routine. The adjustments made to step size are still not enough
to increase the speed of the 5^{th} order routine enough to compete with the
4^{th} order routine.

The error in the energy rises following a pattern similar to the oscillations in the radial direction. This growth of this error in the energy occurs most rapidly when the star is decreasing its radius from the center of the galaxy. Then the error rises at a much slower pace when the radius is increasing. This creates a stepping function in the energy's error.

I would like to take this opportunity to thank Dr. Gonthier for allowing me the opportunity to take part in this program. I would also like to thank the National Science Foundation and the NSF-REU for providing the funding. Hope College and the Laboratory for High Energy Astrophysics at Goddard Space Flight Center for the use of computers and office space.

bj277102@hope.edu