Joel Berrier's pictureJoel Berrier

Hope College

Research Professor - Dr. Peter Gonthier

Supported by the NSF-REU and NSF


Computer Simulation of Galactic Pulsar Distribution

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 5th order Runge-Kutta-Fehlberg method to replace the previously used 4th 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 109 years. The 4th order code also had a crude adaptive step size to help adjust itself to the current error in the energy.

This 5th 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 4th order and 5th order calculation of the relevant quantities. The relevant quantity in our program is the error in the energy. We used the difference between these 4th and 5th 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 5th order Runge-Kutta-Fehlberg routine have been less then what we had desired. The 5th order Runge-Kutta code produces a far larger number of bad events then in the previously used 4th 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 4th order routine also currently runs approximately 2.5 times faster than the recent version of the 5th order code.

The adaptive step size is controlled by the difference between the fractional error in the energy as calculated with both the 5th and 4th order methods. This quantity is scaled to our problem by dividing it with the current energy from the 5th 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 5th 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 5th order routine enough to compete with the 4th 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