TABLE OF CONTENTS

Executive Summary

Chaos and Dynamic Systems 

The Physics

The Mathematics
    Summary of Mathematical Foundations
    Construction of the Equations
    Energy Equations

The Program

Conclusion

References       

Appendix A- Driver Program for Runge-Kutta Solver

Executive Summary:

    We are team 72 from Sandia Preparatory School. The goal of our project has been to model the behavior of the Levitron™, a small magnetic top, marketed as a child’s toy. The Levitron™ is capable of hovering as it spins above a larger magnetic base. The interaction of the magnetic fields of the top and the base, as well as the top’s rotational motion, allow for it to sustain levitation for up to two minutes. We chose this seemingly simple toy as the subject for our project because when one delves into the mathematics and physics behind its behavior, one finds that it actually exhibits chaotic properties. It has also been discovered that subatomic particles trapped in magnetic fields demonstrate behavior similar to that of the Levitron™. We chose to simulate its behavior because our team had a mutual interest in chaos theory and the Levitron™ would give us an introduction to some of the basics behind this newly emerging field such as bifurcation and strange attractors. Chaos theory has become a significant field because it can be applied to virtually any dynamic system including weather patterns or the stock market. In an increasingly complicated world, chaos theory has redefined this complexity and has the potential to become one of the most significant foundations for future scientific studies.
     As most people who have played with magnets before know, it seems to be impossible for one magnet to levitate above another. The Levitron™ accomplishes this feat by virtue of its rotational motion. The rotation of the top creates a gyroscopic motion that keeps it floating as long as there is enough rotational momentum.
    We are modeling its behavior mathematically through Hamilton’s equations of motion, tracking the changes that occur over a given time step. As explained in the body of this paper, we solve 12 differential equations simultaneously using a method known as Runge-Kutta to accommodate the time intervals between steps. Our code outputs sets of values for a single variable. By graphing these variables with respect to time, we have been able to study changes along the ‘x’ and ‘y’ axes and determined that they act as a ‘strange attractor’. We have also modeled the effects of momentum and are currently in the process of examining the chaotic behavior of the angles of rotation and how they relate to each other. The top’s behavior has proved to be an interesting subject that has potential significance for future research on chaos.


Figure 1: Demonstrating the Levitron
in suspension above the magnetic base.

Back to the Table of Contents

Chaos and Dynamic Systems:

    It is important to discuss dynamic systems before delving into the chaotic behavior that many demonstrate. A dynamic system is any system that changes over time. This change can have a distinct pattern or can be completely random. Virtually anything around us that involves change can be classified as a dynamic system. Examples can include everything from weather patterns to traffic patterns, or even the currents in the ocean. Chaotic and random systems emerge mathematically because of their extreme sensitivity to initial conditions. This separates them from more ordered systems; a common example being a Newtonian physical problem in which the answers do not change every time you solve the equations. We can study chaotic systems through graphs that represent the solutions to the system as the initial conditions evolve. Such graphs can include plots of average traffic speed with respect to time, amount of moisture with respect to time, or wind speed with respect to time. Creating such plots allows the system to be studied   graphically and lets us view its patterns. Even such influential and everyday things as the stock market demonstrate chaotic behavior.
    A chaotic system is a type of dynamic system that can be predicted by finding patterns in the solutions to its equations, although a specific solution never repeats. In order to study these systems, we solve several equations over time given certain initial conditions and then plot the results on graphs. The chaotic patterns emerge on these graphs after solving the equations with multiple sets of initial conditions. Changing the conditions helps in the study of the behavior of these systems by allowing us to perceive how each condition affects the system as a whole. The Levitron™ is a good example of this. Factors such as spin rate, height and position along the ‘x’ and ‘y’ axes greatly effect both its stability and behavior.
    There are different types of chaos. “Strange attractors” are one type. When looking at a graph of a chaotic system, for example: precession with respect to time, it will appear as though the graph is trying to get to a certain point but never quite makes it. Strange attractors often result in graphs that look like spirals but without perfect geometric qualities. Bifurcation is another type of chaos that can be seen in the Levitron™. Bifurcation is essentially the multiplication of frequencies, which results in a branching pattern that perpetually increases in complexity. Bifurcation is evident in the Levitron™ upon examining the frequency of the ‘nodding’ or ‘wobble’ of the top itself. This is usually referred to as nutation. The nutation may initially begin with a simple frequency and then appear to become completely random. In fact, additional frequencies are emerging and basically stacking upon one another, making the motion increasingly complex as the top precesses about the spin axis. If more than one nutation frequency is present in our model, it starts to divide exponentially becoming four frequencies, then eight, and so on.

Back to the Table of Contents

The Physics:

    The Levitron™ is a spinning magnetic top that floats above a square base with a ring magnet inside. Certain very precise conditions must be met in order for the system to be stable enough to maintain prolonged levitation. The magnetic field of the base has to be aligned parallel with gravity or else the top will be forced to one side or the other. The rate of spin of the top can neither be too slow nor too fast. If it is spun too slowly, the torque caused from gravity cannot be overcome and the top collapses. If it is spun too quickly, the top maintains an axis of rotation which is too close to the vertical (in other words the tilt is not great enough) and stability cannot be achieved. This tilt is known as precession and is a crucial factor in our model. Without it, the top is prone to collapse. Once stability is achieved, the top will float for a maximum of two minutes because air resistance slows its rotation. The top can sustain levitation for much longer periods in a vacuum. Some vacuum experiments have been carried out, demonstrating that the top will spin for over 30 minutes. The top eventually collapses because the multiplying nutation frequencies increase the angle of precession leading to instability (this will be explained further below). The only way to counteract this is to provide an outside energy source that constantly re-aligns the top to the center of the field. One such device is marketed as the ‘Perpetuator.’ This device keeps the top rotating at the ideal spin rate of about 30 revolutions a second  and prevents the precession angle from becoming too large by the utilization of another rotating magnet. The ‘Perpetuator’ can maintain stable levitation for approximately four hours. The eventual collapse is usually because of changes in the room temperature, which the Levitron™ is sensitive to (this will further explored below as well).
    The height above the base in which the top is stable is very small: about 1.25 to 1.75 inches . Meeting these conditions is accomplished by, first, spinning the top as close to the center of the base as possible so that it can align itself in the center of the magnetic field. Secondly, the lifter plate is brought up slowly until it is at the proper height, at which point the top “jumps” off the plate and stabilizes itself. The weight of the top also has to be very precise in order to balance out the force of gravity and magnetic repulsion. Several thin, metallic washers are used to change the weight of the top, which must be adjusted to within a tenth of a gram. Slight temperature changes to the magnetic base (plus or minus six degrees Fahrenheit) affect the magnet field enough that the weight must be adjusted, thus excessive handling or even changes in room temperature can effect its behavior .
You have seen that the stability of the top is dependent on several factors, however we have not yet examined why these conditions are so sensitive. From Earnshaw’s theorem, we discover that a stationary magnet cannot levitate above another stationary magnet. One property of all bodies is that they prefer to exist in the lowest possible energy state. For stability to be achieved there must be a suitable minimum energy state.


Figure 2.1: Proving that a minimum must be present to achieve
stable levitation. ‘X’=Change in Position, ‘Y’=Energy.
    The first graph in figure 2.1 does not have a true minimum but rather only an inflection point. Mathematically this is because the derivative is zero at that point. Physically, this is because there is literally only one point where the top would be stable. Because this point is infinitesimally small, it is impossible to achieve stability without perfect conditions, which do not exist in reality. However, when the top is spun, it creates a minimum as shown in the second graph, allowing stable levitation.
This minimum is created because of the stabilization caused by gyroscopic motion, which also prevents a normal spinning top from falling over. Since the top cannot be aligned perfectly parallel to gravity, gravity will exert a force on it, tilting it away from the vertical.

Figure 2.2: Precession
    This torque is defined as T = R x F. The cross products of the radius and force (gravity) vectors result in a torque vector that points out from the diagram. Thus, instead of falling over, the top rotates, or precesses, around an axis parallel with gravity. This precession demonstrates chaotic behavior and is one of the factors that we are studying. If the top were perfectly vertical, gravity would not exert a torque on it since the vector would have to be at least partially perpendicular to the axis of rotation. This lack of precession causes the top to start losing stability.

Figure 2.3: Instability is caused when
the top is aligned with the ‘z’ axis
Understanding the basics behind the precession of the top itself (figure 2.3, above) and how its motion lends itself to stability is crucial. Not only is it in this simple precessional model that we discover chaos, but we are also introduced to the primary physical variables that compose the mathematical equations. Once developed, these equations unfold to demonstrate the complexity of what may initially seem to be a simple toy.

Back to the Table of Contents

The Mathematics:
Section 3.1 “A Summary of the Mathematical Foundations”:

    The mathematics involved in the examination of the conditions that govern the stability and motion of the Levitron™ can be described through a set of equations that we shall attempt to describe in the following section. Beginning with certain initial conditions and assumptions, we can mathematically create a stable environment for the spinning top and proceed to solve for the physical characteristics of the system. Essentially, we begin with conditions for height, strength of the dipole, and other such attributes that we know will be conducive to allowing the top to levitate. This allows us to avoid having to 're-discover' exactly where to begin since the region of stability for the top is so narrow.
    Beginning with these conditions, we can move on to examine the equations which stem from the standard equations for the kinetic and potential energy of a spinning body. Using these equations and a method known as nondimensionalization, we construct what is commonly called the Hamiltonian equation. Nondimensionalization allows us to simplify the equations by eliminating certain physical characteristics of the system and treating the top as a single dipole without mass (this will be further explained below). The nondimensionalized form of the Hamiltonian can then, in turn, be used to establish a set of equations that solve for the necessary variables which will be described in more depth. This set of motion equations is analyzed over a given time interval to give values for the variables. The initial conditions are then altered and the process is repeated. The changes to the solutions can be graphed so that we can analyze them visually. It is through this analysis that the chaotic behavior of the Levitron™ becomes evident.

Back to the Table of Contents

Section 3.2 "Construction of the Equations":
    The first important mathematical concept in our problem is the understanding of the two vectors that represent the position and momentum of the system. Position (the ‘q’ vector) contains six variables as shown below.

    The first three variables represent physical coordinates of the spinning top; the last three need further explanation.  is the actual angle of precession, or the angle that the top assumes away from the vertical axis (it could also be referred to as tilt).  represents the arc that the top covers as it rotates and  represents the motion of the top itself around its individual axis (see figure 3.1).

Figure 3.1: Illustrating the Angles associated
with the motion of the top.
Similarly, the ‘p’ vector contains values that represent the same characteristics with respect to time.
    The position coordinates are rather straightforward. We end up with values which are representative of the velocity (position with respect to time). Unlike these position coordinates where the momentum along a given axis is not dependent on the other axes, the angular components (p4, p5, and p6) are more complicated. The factors that describe angular momentum change in more than one direction. This interdependency complicates the ‘p’ vector as you can see above.
    The relation between these vectors should be clear by this point. The ‘q’ vector represents the physical attributes of the system and the ‘p’ vector represents these attributes with respect to time. We identify these components according to their respective vectors and their positions. For instance, the physical coordinate, ‘x’, is q1, ‘y’ is q2, etc. This means that p1 is ‘x’ taken with respect to time and so on. These two vectors contain all the variables we need to describe the Levitron™. From here, we can move on and apply these attributes mathematically.

Back to the Table of Contents

Section 3.3 "Energy Equations”:
    Once we have established the two primary vectors in our problem, we can begin to construct the essential equations that will allow us to evaluate these variables.
    We begin with the equations for kinetic and potential energy for a rotating body, which we are familiar with from basic physics. We know that kinetic energy is equal to 1/2 the mass times the velocity squared (1/2mv^2) in the simplest mechanical context. In this case, we expand this to contain the necessary attributes of the spinning top:

    You can see that this equation can be dissected to return to the basic form that was just mentioned. Since we are taking the position with respect to time (X, Y, etc.), we can substitute this for velocity, as basic calculus dictates. The latter portion of this equation in which we see that
represent the transverse and polar moments of inertia respectively. The transverse moment of inertia is perpendicular to the main spin access whereas the polar moment of inertia is along the axis itself. The potential equation is as follows:
    We mentioned earlier that we would be employing a method known as nondimensionalization for the simplification of our equations. The equations for energy above do not demonstrate this since they contain the variable ‘m’, or mass. By treating the base as a ring-dipole and the spinning top as a single dipole without mass, we can eliminate some of the variables in the equations and simplify the solution. We take the variables that require dimensions and represent them as constants. In our problem the constants ‘a’, ‘c’, and ‘’ contain the dimensional variables, including mass. The values chosen for the constants are based on scientific observations of the Levitron™’s behavior. Nondimensionalization is one way to arrange the essential elements while eliminating a clutter of unnecessary variables and representing them as more efficient constants. This not to say that the results are any less accurate, we are merely utilizing this particular numerical method for simplicity’s sake.
 From these two equations, we can compile what is known as the Hamiltonian equation. This will be used in Hamilton’s equations of motion that will be presented below. By combining the kinetic and potential descriptions of the system, we can effectively create a single equation that allows us to evaluate each of the previously mentioned vectors with respect to each other. The Hamiltonian equation in its pure form is shown below:
    Once more, the kinetic portion is evident in the first three terms (p1^2, p2^2, p3^2). Again, recall that the first three ‘p’ components squared are equal to velocity since they are the physical coordinates, or position, of the system taken with respect to time. The second two terms are the transverse moment of inertia while the p26/c is the polar moment of inertia. In the second half of the equation, you will notice the aforementioned constant ‘mu’. This represents the strength of the dipole and thus is the magnetic potential energy. This is dependent on the tilt of the top as well as the other angular attributes. The final term, q3, is essentially the potential energy after nondimensionalization. We are familiar with potential energy as mgh, but if we remove mass and gravity (since gravity is negligible for a single dipole without mass in a magnetic field) we are left with the height, or q3. If you look back to the presentation of the ‘q’ vector (equation 1), you will see that q3 is ‘z’, or height on a physical axis system. Although the purpose of this equation may still be rather unclear initially, it is applied in a concise and streamlined manner in Hamilton's equations of motion. Using the two equations of motion, shown below, we can evaluate either vector based on the partial derivative of the Hamiltonian equation above with respect to the other vector. The following two equations represent what is finally derived from the compilation of the equations above:
    Thus position can be evaluated according to the partial derivative of the Hamiltonian equation with respect to the ‘p’ vector and vice-versa. We actually only evaluate one variable at a time. For instance, q1 is equal to the partial derivative of the Hamiltonian with respect to p1 and so on.
    This leaves us with twelve differential equations, one for each component in each of the two vectors. We can solve these twelve equations simultaneously and from this we can successfully model the Levitron™. If we begin with conditions for one vector, we can evaluate the other to end up with a new set of values. Continuing to solve over a given time step allows us to continually track and analyze the characteristics of the system as the top precesses through its rotation. The following section will attempt to explain, in depth, how we use a method known as Runge-Kutta to solve our motion equations repeatedly over the given time step. Acquiring new sets of values for the vectors each time will provide us with data that we can plot and analyze graphically. As explained above in the section titled ‘Chaos and Dynamic Systems’, this is where the chaotic behavior of the Levitron™ becomes evident and it is here that our project culminates.

Back to the Table of Contents

The Program:

     In order to solve the twelve equations presented above, we will be utilizing the Runge-Kutta method. Runge-Kutta works by transforming the differential equation in question to its slope form. Both sides of the equation are multiplied by the change in ‘X’ or, in our case, the change in ‘T’ (time). This allows us to solve for the current variable by plugging in a step size for the change in ‘T’ and a previous value for the variable we are solving for. This process is repeated to obtain sets of values for the variable.
    This is the simplest form of the Runge-Kutta method and is known as a ‘First-Order Runge-Kutta Method’ or ‘Euler’s Method’. This particular method has many limitations including a tendency to become unstable and a relatively high error. To counter this, Runge-Kutta methods have been developed that use more than one previous value to predict the behavior of the equation in the future. This approach is called a ‘Higher-Order Runge-Kutta Method’. This has the benefit of being relatively stable and fairly accurate. The problem is that for each previous value that we consult to get the current one, we add computation time. With the relatively small time step we are using for our problem, this computation time adds up quickly; this is why we need the super computers.
    Runge-Kutta relies heavily on initial conditions. For the first time step we must define each variable. As we have said before, the goal of our project was not to study when the Levitron™ is stable but to study its chaotic behavior. We want initial conditions that provide for stable and continuos levitation. Several papers have been written about the Levitron™ already, focusing mainly on stability. We can use their results for the initial conditions and constants to create an environment suitable for stable levitation. This allows us to build on previous results and study the chaotic behavior -- an area left relatively unexplored. The keys to stable levitation are the height of the top and the momentum of the rotation. The Levitron™ is only stable at one height and with a high enough momentum to counteract the torque from the base magnet. If we do not meet these conditions, the system will collapse. Another interesting factor is that if the momentum is too high, the system will also destabilize; since the user controls this entirely, it is difficult to achieve just the proper speed and twist when releasing the top. This only reiterates how difficult it is to create stable levitation if the physical parameters are not met. The height of the top, q3, is initialized at 1.72 inches. The tilt, q4, is small (usually initialized between 0.001 and 0.015 radians) . The momentum of the top itself, p6, is set between 5.0 and 9.0 radians per second . This momentum is constant. The last value we must initialize is the momentum around the vertical axis, p5, which is usually between 0.5 and 3.5 radians per second . We fixed the physical constants at a=0.089, c=0.139, and =8.20; which are consistent with observed values for the Levitron™ . The final facet to be decided, after the initial conditions and constant values have been dealt with, is the time step. With the nondimensionalization of the equation, the consistent time scale is equal to   which is the time it would take the top to fall a distance of R/2 (for our model R=34.7mm so time is equal to 59.5ms.) One thousand units of this time scale are equal to approximately one minute.
    One final note about our model: We do not take into account the effects of air resistance. The only thing that air resistance does is to slow down the rate of spin leading to collapse of the system. We can reasonably assume that air resistance does not change the chaotic behavior of the top in any great degree. Since our goal was to study the chaotic behavior, we did not see the necessity of adding unneeded complexity.
    We used several different Runge-Kutta routines, as well as different computers. The first method was in ‘Waterloo Maple V’. We programmed the equations into Maple V in order to solve small problems on a PC that could be used to verify the results we receive from our actual C code. The Maple V code functioned well but we were unable to solve anything but the simplest problems due to the PC’s inability to process a low enough error tolerance. The second option we used was from the book, Numerical Recipes in C. We have programmed a driver that utilizes the 5th order Runge-Kutta solver available in this book. The results obtained from the program are compared to the simple problem we solved using Maple V, and several test cases were compared to previously reported results.
    We are graphing the data using ‘Dataplot’ from Fortner software. We have included three of these plots. The first two demonstrate how a low initial momentum causes instability. Figure 4.1 shows a graph of tilt versus time units. Remember that our time unit is 59.5ms so the graph represents about half a second of real time. Initially the precession of the top is within a very small range (between 0.001 to 0.015 radians) but the spin rate is too low to counter the torque from gravity, causing the tilt angle to increase at the 8th time unit. The increased tilt angle is too high for the Levitron™ to continue levitating and the top collapses. Figure 4.2 shows a graph of height versus time units and you can see how the collapse of the top matches up with the increase in the tilt angle. The last graph represents part of the top’s chaotic behavior. Figure 4.3 is a graph of the changes along the ‘x’ and ‘y’ axes. The chaotic behavior that this position graph demonstrates shows one of the types mentioned in the ‘Dynamics section,’ a strange attractor. If you look carefully you can see the spiral-like change in the solutions. The graph represents 165 time units (approximately 10 seconds). We initialized the momentum of the top, p6, to 8.7 radians a second; this is very close to the maximum value under stable levitation. We are currently studying the chaotic behavior of the angles and how they relate to time as well as to each other. These relationships are where the truly interesting chaos emerges.


Figure 4.1: This graph shows tilt versus time

Figure 4.2: This shows height versus time.

 Figure 4.3: This shows the changes in the ‘x’ and ‘y’ position.

Back to the Table of Contents

Conclusions:
     We have successfully modeled the Levitron™ for both unstable and stable configurations. We verified are results with previously reported results on the stability of the Levitron™ and by comparing the results to simpler solutions completed using Maple V on a PC.
    Our most original accomplishment this year has been the study of the chaotic behavior of the Levitron™. As we said earlier in the paper, studies have been conducted on the Levitron™ in the past but have focused almost exclusively on stability. We have primarily looked at the chaotic behavior, building upon most previous work.
    What remains to be done includes further study on the chaotic behavior between the angles of rotation. We are currently attempting to model the Levitron™ during sustained levitation in order to graph the changes in the angles. The chaotic behavior emerges in these systems when they are taken with respect to time and with respect to each other.
    We have learned a great deal over the course of the project, including how to work with non-linear mathematics. In addition, we have gained a lot of practice in the area of computing with initial value problems. We have come to understand the basics of chaotic theory as well as the value in studying this field.
    We would like to thank everyone at Los Alamos National Labs for providing us the resources to complete these projects and for their support over the course of our high school careers. It has been a very positive program.
    We owe the most gratitude to our advisor of four years, Dr. Kevin Malloy, who has given up countless hours to help us understand the math and science behind our projects. We cannot forget his colleagues at the Center for High-Technology Materials at the University of New Mexico who have provided technical assistance as well as access to software packages. We extend a special thanks to Dr. Fred Gelbard from Sandia National Laboratories who helped us debug our program and who was there to answer questions when our advisor was not available. Last but certainly not least, we wish to thank our sponsor, Mr. Neil McBeth who originally introduced us to the program four years ago and who has always been there for us.

Back to the Table of Contents

References:

Berry, Michael V. The Physics of the Levitron. (Bibliographic data unknown. From the
        pamphlet that is included with the Levitron.)

Chua, Leon; Madan, Rabinder N.. Sights and Sounds of Chaos. New York: IEEE Service
        Center, 1988. (Originally published in IEEE Circuits and Devices Magazine).

Dullin, Holger R.; Easton, Robert W.. Stability of the Levitron. Boulder, Colorado:
        University of Colorado, April 23, 1998.

Gans, Roger; Jones, Thomas B.; Washizu, Masao. Dynamics of the Levitron. UK: IOP
        Publishing Ltd., 1998.

Press, William; Flannery, Brian; Teukolsky, Saul; Vetterling, William. Numerical Recipes
        in C: The Art of Scientific Computing. Cambridge: Cambridge University Press, 1988.

 Simon, Martin D.; Heflinger, Lee O.; Ridgeway, S.L..Spin Stablized Magnetic Levitation.
    California: UCLA Dept. of Physics (published in the American Journal of Physics), unknown pub. date.

 As well as the following software packages:
  Waterloo Maple V
  Fortner Software’s “DataPlot”

Back to the Table of Contents

Appendix A

Driver Program for Runge-Kutta Solver

#include <stdio.h>
#include "nr.h"
#include "nrutil.h"
#include <math.h>

#define N 12 /*The number of simultaneous equations to solve*/
#define ITER 100
#define a 0.089
#define c 0.139
#define M 8.1880
#define MAXSTP 10000
#define TINY 1.0e-30

int kmax=0,kount=0;  /* defining declaration */
float *xp=0,**yp=0,dxsav=0;  /* defining declaration */

void derivs(x,y,dydx) /*Function to solve the right hand side of the equations                               used by odeint. DyDx is an array that holds the answers for each variable
represented the y array.*/
float x,y[],dydx[];
{
 dydx[1] = y[7];
 dydx[2] = y[8];
 dydx[3] = y[9];
 dydx[4] = (y[10]/a);
 dydx[5] =0
(y[11]-(y[12]*cos(y[4])))/(pow(sin(y[4]),2.0));
 dydx[6] =
(y[12]/c)-(((y[11]-(y[12]*cos(y[4])))*((cos(y[4])/sin(y[4]))*(1/sin(y[4]))))/a);
 dydx[7] =
(1/(4*pow(((1.0+pow(y[3],2.0))),(9.0/2.0))))*(3*M*(-y[1]*cos(y[4])*(3-24*pow(y[3],2.0)
+8*pow(y[3],4.0))+y[3]*cos(y[5])*sin(y[4])*(-3-pow(y[3],2.0)+2*pow(y[3],4.0))));
 dydx[8]=
(1/(4*pow((1.0+(pow(y[3],2.0))),(4.5))))*(3*M*(-y[2]*cos(y[4])*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0))
+y[3]*sin(y[5])*sin(y[4])*(-3-pow(y[3],2.0)+2*pow(y[3],4.0))));
 dydx[9]=
(1/(8*pow((1.0+(pow(y[3],2.0))),(5.5))))*(-3*M*y[3]*cos(y[4])
*(-12-75*pow(y[2],2.0)-16*pow(y[3],2.0)+200*pow(y[2],2.0)*pow(y[3],2.0)
+4*pow(y[3],4.0)-40*pow(y[2],2.0)*pow(y[3],4.0)+8*pow(y[3],6.0)-5*pow(y[1],2.0)
*(15-40*pow(y[3],2.0)+8*pow(y[3],4.0))+(2*(1.0+pow(y[3],2.0))
*(2*pow((1.0+pow(y[3],2.0)),(.5))+8*pow(y[3],2.0)*pow((1.0+pow(y[3],2.0)),(.5))
+12*pow(y[3],4.0)*pow((1.0+pow(y[3],2.0)),(.5))+8*pow(y[3],6.0)*
pow((1.0+pow(y[3],2.0)),(0.5))+2*pow(y[3],8.0)*pow((1.0+pow(y[3],2.0)),(0.5))
-3*M*y[1]*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0))*cos(y[5])*sin(y[4])
-9*M*y[2]*sin(y[4])*cos(y[5])+2*pow(y[3],8.0)*pow((1.0+pow(y[3],2.0)),(0.5))
-3*M*y[1]*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0))*cos(y[5])*sin(y[4])
-9*M*y[2]*sin(y[4])*sin(y[5])+72*M*y[2]*pow(y[3],2.0)*sin(y[4])*sin(y[5])
-24*M*y[2]*pow(y[3],4.0)*sin(y[4])*sin(y[5])))));
 dydx[10]=
(1/2)*(((2*y[12]*(y[11]-y[12]*cos(y[4]))/a*sin(y[4]))*((2*pow((y[11]
-y[12]*cos(y[4]),2.0)*cos(y[4]))/a*pow(sin(y[4]),2.0)*sin(y[4]))
-(1/(4*pow((1.0+(pow(y[3],2.0))),(4.5))))*(M*((-4-9*pow(y[2],2.0)+72*pow(y[2],2.0)
*pow(y[3],4.0)+8*pow(y[3],6.0)-3*pow(y[1],2.0)*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0)))
*sin(y[4])-6*y[3]*(-3-pow(y[3],2.0)+2*pow(y[3],4.0)*cos(y[4])*(y[1]*cos(y[5])
+y[1]*sin(y[5]))))))));
 dydx[11]=
(3*M*y[3]*(-3+2*pow(y[3],2.0))*sin(y[4]))*(y[2]*cos(y[5])-y[1]*sin(y[5]))/
((4*pow((1.0+(pow(y[3],2.0))),(3.5))));
 dydx[12]=0;
}

void odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqc) /*Driver for Runge-Kutta solver*/

float ystart[],x1,x2,eps,h1,hmin;
int nvar,*nok,*nbad;
void (*derivs)(); /* ANSI: void (*derivs)(float,float *,float *); */
void (*rkqc)();  /* ANSI: void (*rkqc)(float *,float *,int,float *,float,
    float,float *,float *,float *,void (*)()); */
{
 int nstp,i;
 float xsav,x,hnext,hdid,h;
 float *yscal,*y,*dydx,*vector();
 void nrerror(),free_vector();
 
 yscal=vector(1,nvar); /*Define arrays so that they start at 1. This is       necessary because the routines were ported from
     Fortran.*/
 y=vector(1,nvar);
 dydx=vector(1,nvar);
 x=x1;
 h=(x2 > x1) ? fabs(h1) : -fabs(h1);
 *nok = (*nbad) = kount = 0;
 for (i=1;i<=nvar;i++) y[i]=ystart[i];
 if (kmax > 0) xsav=x-dxsav*2.0; /*Guarantee storage of first step*/
 for (nstp=1;nstp<=MAXSTP;nstp++) { /*Take MAXSTP at the most*/
  (*derivs)(x,y,dydx);
  for (i=1;i<=nvar;i++)  /*Scale accuracy*/
   yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
  if (kmax > 0) {
if (fabs(x-xsav) > fabs(dxsav)) {  /*Store intermediate values*/
    if (kount < kmax-1) {
     xp[++kount]=x;
     for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
     xsav=x;
    }
   }
  }
  if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x; /*Cut down step size?*/
  (*rkqc)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,derivs);
  if (hdid == h) ++(*nok); else ++(*nbad);
  if ((x-x2)*(x2-x1) >= 0.0) {   /*Are we done?*/
   for (i=1;i<=nvar;i++) ystart[i]=y[i];
   if (kmax) {
    xp[++kount]=x;
    for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
   }
   free_vector(dydx,1,nvar);
   free_vector(y,1,nvar);
   free_vector(yscal,1,nvar);
   return;     /*Normal exit*/
  }
  if (fabs(hnext) <= hmin) nrerror("Step size too small in ODEINT");
  h=hnext;
 }
 nrerror("Too many steps in routine ODEINT");
}

main()
{
 int i,nbad,nok,j;
 float eps,h1,hmin,x1,x2,*ystart;

 ystart=vector(1,N);

 x1=0.0;  /*Initial time step*/
 x2=0.002;

 ystart[1]=0.0;  /*Set the initial values for the variables*/
 ystart[2]=0.0;
 ystart[3]=1.72;
 ystart[4]=100.0;
 ystart[5]=10.0;
 ystart[6]=0.0;
 ystart[7]=0.0;
 ystart[8]=0.0;
 ystart[9]=0.0;
 ystart[10]=0.0;
 ystart[11]=3.49;
 ystart[12]=4.8;

 eps=1.0e-4;    /*Error tolerance*/

 h1=1.0;
 hmin=0.0;
 
 for(j=1;j<=ITER;j++) {
 
  x1=x1+0.002; /*Increase time step*/
  x2=x2+0.002;
 
odeint(ystart,N,x1,x2,eps,h1,hmin,&nok,&nbad,derivs,rkqc); /*Call Runge-Kutta
Driver*/
 
  for(i=1;i<=N;i++) {
   printf("%f ",ystart[i]);
  }
 
  printf("\n");

 }
 free_vector(ystart,1,N);
}

Back to the Table of Contents