Verlet algorithm: Why am I getting this output?

In summary, my code produces no acceleration in the ##y##-direction (aside from for ##t_n = 0##,) and that the planet just flies off to infinity without finding a stable orbit.
  • #1
H Smith 94
Gold Member
55
1
Hi! I am currently trying to write a code in C to simulate the orbit of planets around the Sun in the solar system.

I am using the velocity Verlet approach and finding that my code produces no acceleration in the ##y##-direction (aside from for ##t_n = 0##,) and that the planet just flies off to infinity without finding a stable orbit.

C:
   /*   INITIALIZATION */

   int N = 366;       // Maximum number of steps taken
   if(N >= MAX_STEPS)   {   // Error handling
     printf("N = %3i is greater than the alloted array space: %4i. \nPlease use smaller N or increase size of arrays.", N, MAX_STEPS);
     exit(0);
   }

   double dt = 1.0 / 365.4;   //   Integration step

   /*Initial conditions*/
   t0     = 0.0;
   t[0]    = t0;
   t[-1]    = t0 - dt;

   r[0]    = 1.0;
   th[0]   = 45*deg;

   x[0]    = sqrt(pow(r[0],2)/(1+pow(tan(th[0]),2)));         // Setting initial position of planet
   y[0]    = sqrt(pow(r[0],2)*pow(tan(th[0]),2)/(1+pow(tan(th[0]),2)));

   double v0;
//   v0   = 2.0*M_PI*Au/yr;
   v0    = sqrt(G*M/r[0]);
   vx0   = v0*cos(th[0]);
   vy0   = v0*sin(th[0]);

//   ax[0] = - (G*M/pow(r[0],2))*cos(th[0]);
//   ay[0] = - (G*M/pow(r[0],2))*sin(th[0]);

   x[-1] = x[0] - vx0*dt;           // Defining x for pre-leap
   y[-1] = y[0] - vy0*dt;           // Defining y for pre-leap

   printf("\tr(-1) = (%f, %f)\n",
       x[-1], y[-1]);

   /********************/
   /* VERLET ALGORITHM */
   /*__________________*/

   int n;                   //   Iteration variable
   FILE *f = fopen(filename, "w");       //   Opening output file in "write" mode
   for(n = 0; n <= N; n++)     {

     r[n] = sqrt(pow(x[n],2) + pow(x[n],2));

     ax[n] = - G*M*cos(th[n])/pow(r[n],2);
     ay[n] = - G*M*sin(th[n])/pow(r[n],2);

     x[n+1] = 2*x[n] - x[n-1] + ax[n]*pow(dt,2);         //   Calls function to step x from x[n] -> x[n+1]
     y[n+1] = 2*y[n] - y[n-1] + ay[n]*pow(dt,2);

     fprintf(f,"%3.4f\t%3.4f\t%3.4f\n",
       x[n], y[n], t[n]);         // Prints plot variables (x,v,t) to file
     printDiagnostics(x, y, ax, ay, t, n);   // Calls subroutine to print diagnostics to Console

     t[n+1] = t[n] + dt;           //   Defines time at iteration n; discrete time Tn = n*dt
   }
   fclose(f);                 // Closing output file

   return 0;                 // To satisfy integer requirement of function
}

This code produces the following output graph:

verletPlot_broken-3D.png

Fig. 1
Output from broken simulation showing no stable planetary orbit.
Note that the Sun is placed at (0,0) and we are examining the Earth's orbit.
I've also provided the full output file from the software as an attachment; the first few lines of diagnostics output are as follows:

Code:
         r(-1) = (0.601989, 0.792215)
n = 0 | r(0) = (0.591813, 0.806075)   -> 0.836950    | a(0) = (33.353809, -45.429385)   t(0) = 0.000000
n = 1 | r(1) = (0.581886, 0.819596)   -> 0.822911    | a(1) = (-58.298016, -0.000000)   t(1) = 0.002737
n = 2 | r(2) = (0.571523, 0.833117)   -> 0.808255    | a(2) = (-60.431386, -0.000000)   t(2) = 0.005473
n = 3 | r(3) = (0.560707, 0.846637)   -> 0.792959    | a(3) = (-62.785279, -0.000000)   t(3) = 0.008210
  ##\vdots##                      ##\vdots##                ##\vdots##                    ##\vdots##                ##\vdots##              ##\vdots##

Is anyone able to see where my code is going wrong? Clearly these graphs do not show a stable orbit!

I am using the MinGW GCC compiler through Eclipse Mars on Windows 8.

It's probably worth noting that this is a homework question, so I'm looking mainly for hints and advice. Thanks in advance!
 
Last edited:
Technology news on Phys.org
  • #2
Have you tried changing your timestep?
If it's too large you'll lose a lot of accuracy.

If I recall correctly I used around 90000 seconds (25 hrs) when I did a similar thing for Halley's comet.
It should be noted however that the period of the comet is much larger, about 75 years.

The second thing I noticed is that your algorithm doesn't look like the velocity verlet algorithm.
When running into problems like this I would always start from the most explicit version of the algorithm.
You can try to simplify it, although I doubt it can get any simpler.
 
  • #3
Just out of curiosity - why are you using negative index values like x[-1]? Did you increment x as a pointer and then want to get the value x[0]? It is not illegal syntactically but you can possibly get undefined behavior and/or weird answers if the code doesn't segfault (raise a memory access exception).

Example:
Code:
 x[-1] = x[0] - vx0*dt;           // Defining x for pre-leap
   y[-1] = y[0] - vy0*dt;           // Defining y for pre-leap
 
  • #4
The immediate error is here:
C:
     ax[n]=- G*M*cos(th[n])/pow(r[n],2);
     ay[n]=- G*M*sin(th[n])/pow(r[n],2);

You are never setting th[n]. Given the behavior, the array is filled with zeroes. Since ##\sin(0)=0##, you are getting no acceleration in the y direction. You should either compute each th[n], or better yet, use a different algorithm. You should instead be using ##\vec a(t) = -\frac {\mu}{||\vec r(t)||^3} \vec r(t)## calculate your accelerations, where ##\mu## is conceptually the equivalent of ##GM##. In code,
C:
     rsq = x[n]*x[n] + y[n]*y[n];
     r[n] = sqrt(rsq);
     fact = -mu/rsq;
     ax[n] = fact*x[n];
     ay[n] = fact*y[n];
     x[n+1] = 2*x[n] - x[n-1] + ax[n]*dtsq;
     y[n+1] = 2*y[n] - y[n-1] + ay[n]*dtsq;

I'm using some auxiliary variable in the above. rsq is the square of the radius, ##||\vec r(t)||^2 = x^2 + y^2##, fact is ##-\frac{\mu}{||\vec r(t)||^3}##, mu is ##\mu## (see above and below), and dtsq is the square of the time step, calculated once outside the loop.

Note also that I don't call pow. That's an old habit. In C and C++, it's better to use x*x instead of pow(x,2).

Regarding the use of mu. The gravitational constant ##G## is only known to four places of accuracy. The standard gravitational parameter of the Sun is known to more than ten places.
 
Last edited by a moderator:
  • Like
Likes H Smith 94, Silicon Waffle and jim mcnamara
  • #5
Thank you for your responses!

JorisL said:
Have you tried changing your timestep?
If it's too large you'll lose a lot of accuracy.

If I recall correctly I used around 90000 seconds (25 hrs) when I did a similar thing for Halley's comet.
It should be noted however that the period of the comet is much larger, about 75 years.

I was aware is had to be small but I didn't know it was that small! I've decreased it to about 1/10th of a day and it has indeed improved the accuracy. I am still getting similarly strange output, as follows.

verletPlot_broken-3D1.png

Fig. 2
Output with timestep ##dt = 1/3650##.​
JorisL said:
The second thing I noticed is that your algorithm doesn't look like the velocity verlet algorithm.
When running into problems like this I would always start from the most explicit version of the algorithm.
You can try to simplify it, although I doubt it can get any simpler.

Honestly, I'm quite a newbie at both C programming and computational physics. I attempted to implement the second algorithm listed on that page: [tex]\mathbf{r}(t+\Delta t) = \mathbf{r}(t) + \mathbf{v}(t) \Delta t + \frac{1}{2} \mathbf{a} \Delta t[/tex] [tex]\mathbf{a}(x+\Delta t) = \frac{1}{m} \mathbf{F}(t+ \Delta t)[/tex] where ##\mathbf{F}(t+ \Delta t)## is Newton's law of gravitation; and [tex]\mathbf{v} (t+\Delta t) = \mathbf{v} (t) + \frac{1}{2}(\mathbf{a} (t) + \mathbf{a} (t + \Delta t)\,\Delta t.[/tex]
I interpreted this as defining the timestep of the next iteration as ##t_n+1 = t_n + \Delta t## and then replacing ## t + \Delta t ## with the relevant array element ## n+1 ##. Clearly, this hasn't worked -- what would be a better method to 'translate' these equations into a code?

UPDATE:

I have just realized I was in fact not using the above algorithm -- thank you for pointing this out! In this particular code I was using the Verlet equation from my lecture notes (pg. 11). This also answers the following question.

Sorry, I have re-written this code using a million different methods and they're all starting to blur together at this point!

————————————————————————

jim mcnamara said:
Just out of curiosity - why are you using negative index values like x[-1]? Did you increment x as a pointer and then want to get the value x[0]? It is not illegal syntactically but you can possibly get undefined behavior and/or weird answers if the code doesn't segfault (raise a memory access exception).
 
Last edited:
  • #6
Wow, thank you for your extensive post!

D H said:
The immediate error is here:
C:
     ax[n]=- G*M*cos(th[n])/pow(r[n],2);
     ay[n]=- G*M*sin(th[n])/pow(r[n],2);

You are never setting th[n]. Given the behavior, the array is filled with zeroes. Since ##\sin(0)=0##, you are getting no acceleration in the y direction. You should either compute each th[n], or better yet, use a different algorithm.

Well spotted! I'm a little embarrased I didn't notice that myself. I have now defined th[n] as th[n+1] = atan( y[n+1] / x[n+1]), or [tex] \theta(t+\Delta t) = \mathrm{arctan}\left(\frac{y(t+\Delta t)}{x(t+\Delta t)}\right).[/tex]
UPDATE:
The following step makes this line of code not-needed. No need to use that pesky atan() function.

D H said:
You should instead be using ##\vec a(t) = -\frac {\mu}{||\vec r(t)||^3} \vec r(t)## calculate your accelerations, where ##\mu## is conceptually the equivalent of ##GM##. In code,
C:
     rsq = x[n]*x[n] + y[n]*y[n];
     r[n] = sqrt(rsq);
     fact = -mu/rsq;
     ax[n] = fact*x[n];
     ay[n] = fact*y[n];
     x[n+1] = 2*x[n] - x[n-1] + ax[n]*dtsq;
     y[n+1] = 2*y[n] - y[n-1] + ay[n]*dtsq;

I'm using some auxiliary variable in the above. rsq is the square of the radius, ##||\vec r(t)||^2 = x^2 + y^2##, fact is ##-\frac{\mu}{||\vec r(t)||^3}##, mu is ##\mu## (see above and below), and dtsq is the square of the time step, calculated once outside the loop.
Should I instead be using ax[n] = fact*x[n]/r[n] here to account for the vector components?

D H said:
Note also that I don't call pow. That's an old habit. In C and C++, it's better to use x*x instead of pow(x,2).
I was not aware of this, excellent point! I had figured I should be a little wary of C's built-in functions! :angel:

D H said:
Regarding the use of mu. The gravitational constant ##G## is only known to four places of accuracy. The standard gravitational parameter of the Sun is known to more than ten places.
Thank you, these are all great pointers! I have updated my code to reflect this and will keep this in mind in future.

My code now outputs the following, which I assume will start working with the correct selection of initial conditions?

verletPlot_broken-3D-2.png

Fig. 3
Output with updating th[n].​
 
Last edited:
  • #7
I've added my loop from the MATLAB code I have laying around. (it's almost the same as C, I will assume you get the gist of it)
You'll notice I never used the goniometric functions meaning I circumvented the issue D H noticed with your code.

First let me clarify the velocity-verlet algorithm (VVM).

A naive way to integrate the equations of motion (which is what we are actually doing) is the one you are using.
It's called the Störmer-Verlet algorithm.
Basically you give some initial values and then start iterating using the rule
##x_{n+1}=2x_n- x_{n-1}+ a^{(x)}_n\Delta t##
Here I defined ##a^{(x)}## as the x-component of the acceleration.
The velocity has been effectively eliminated from this method.

The VVM however uses the velocity in a very explicit way.
The algorithm consists of some extra steps.

  1. Calculate ##v^{(x)}(t+\frac{1}{2}\Delta t) = v^{(x)}(t) + a^{(x)}(t) \frac{1}{2}\Delta t##
  2. Calculate ##x(t+\Delta) = x(t) + v^{(x)}(t+\frac{1}{2}\Delta t)##
  3. Find ##a(t+\Delta t)## using the interaction potential, in this case the gravitational potential
  4. Update ##v(t+\Delta t) = v^{(x)}(t+\frac{1}{2}\Delta t) + a^{(x)}(t+\Delta t) \frac{1}{2}\Delta t##
Notice how I neglected to save the first step and plugged it directly into the second step.
This isn't a very good idea, since I calculate the same thing twice each iteration.
For this simple example this doesn't give us much if any trouble. For larger simulations you want to avoid that.

Matlab:
for ii = 2:N
  
   r(ii,1) = r(ii-1,1)+v(ii-1,1)*Dt+Ax/2*Dt^2;
   r(ii,2) = r(ii-1,2)+v(ii-1,2)*Dt+Ay/2*Dt^2; 
   R = sqrt(r(ii,1)^2+r(ii,2)^2);

   %First velocity update
   v(ii,1) = v(ii-1,1)+Ax/2*Dt;
   v(ii,2) = v(ii-1,2)+Ay/2*Dt;

   %Update acceleration
   Ax = (-G*M/R^2)*(r(ii,1)/R);
   Ay = (-G*M/R^2)*(r(ii,2)/R);

   %Second velocity update
   v(ii,1) = v(ii,1)+Ax/2*Dt;
   v(ii,2) = v(ii,2)+Ay/2*Dt;
  
   %Energy
    K(ii) = 1/2*m*(v(ii,1)^2+v(ii,2)^2);
    V(ii) = -G*m*M/R;
end

EDIT: (crosspost)
Your latest figure seems to point towards a division by a very small number.
 
  • Like
Likes H Smith 94
  • #8
JorisL said:
I've added my loop from the MATLAB code I have laying around. (it's almost the same as C, I will assume you get the gist of it)
This is a very useful resource! I'm a little wary about implementing half-steps at this point but it has certainly helped me relate the code to the math, thank you.

JorisL said:
Your latest figure seems to point towards a division by a very small number.
I thought the same thing, too! But I can't seem to find where it would be, unless 1.0 Au counts as a small number for the initial position r[0]. If you can see anything blaringly obvious, my updated code is as follows:
C:
  /********************/
   /* VERLET ALGORITHM */
   /*__________________*/

   int n;                   //   Iteration variable
   FILE *f = fopen(filename, "w");       //   Opening output file in "write" mode
   for(n = 0; n <= N; n++)     {
   
     r[n] = sqrt(pow(x[n],2) + pow(y[n],2));
   
     double dtsq = dt*dt;
     double rsq  = r[n]*r[n];
     double xhat = x[n]/r[n];
     double yhat = y[n]/r[n];

     ax[n] = - (G*M/rsq)*xhat;
     ay[n] = - (G*M/rsq)*yhat;

     x[n+1] = 2*x[n] - x[n-1] + ax[n]*dtsq;         //   Calls function to step x from x[n] -> x[n+1]
     y[n+1] = 2*y[n] - y[n-1] + ay[n]*dtsq;

//     th[n+1] = atan( y[n+1] / x[n+1]);

     fprintf(f,"%3.4f\t%3.4f\t%3.4f\n",
       x[n], y[n], t[n]);         // Prints plot variables (x,v,t) to file
     printDiagnostics(x, y, ax, ay, t, n);   // Calls subroutine to print diagnostics to Console

     t[n+1] = t[n] + dt;           //   Defines time at iteration n; discrete time Tn = n*dt
   }
 
Last edited:
  • #9
Thank you for all your help everyone!

verletPlot_2D_EARTH.png

Fig. 4
Complete orbit of the Earth around the Sun in 1 year.​

I have attached the full code as a .txt file for anyone interested.
 

Attachments

  • VerletSim_C.txt
    7.3 KB · Views: 527

Related to Verlet algorithm: Why am I getting this output?

1. What is the Verlet algorithm and how does it work?

The Verlet algorithm is a numerical method used for integrating equations of motion in classical mechanics. It is a popular algorithm for simulating the motion of particles in molecular dynamics simulations. The algorithm works by calculating the position and velocity of particles at a given time step based on their previous positions and accelerations.

2. Why is the Verlet algorithm used instead of other integration methods?

The Verlet algorithm is often used because it is more accurate and stable compared to other integration methods, such as the Euler method. It also conserves energy, making it a better choice for long-term simulations.

3. What could cause unexpected or incorrect output when using the Verlet algorithm?

There are several reasons why you may be getting unexpected or incorrect output when using the Verlet algorithm. Some common causes include incorrect initial conditions, errors in the force calculations, or incorrect implementation of the algorithm.

4. How can I troubleshoot issues with the output when using the Verlet algorithm?

If you are experiencing issues with the output of the Verlet algorithm, it is important to first check your code for any errors or mistakes. You can also try adjusting the time step or checking the forces acting on your particles. If all else fails, it may be helpful to consult with other researchers or experts in the field.

5. Are there any limitations or drawbacks to using the Verlet algorithm?

Like any numerical method, the Verlet algorithm also has some limitations and drawbacks. One major limitation is that it is only accurate for systems with conservative forces. It may also have difficulty handling highly nonlinear systems. Additionally, the algorithm can be computationally expensive for large systems or long simulations.

Similar threads

  • Programming and Computer Science
Replies
9
Views
1K
  • Programming and Computer Science
Replies
4
Views
1K
  • Programming and Computer Science
Replies
3
Views
1K
  • Programming and Computer Science
Replies
31
Views
2K
  • Programming and Computer Science
Replies
1
Views
692
  • Programming and Computer Science
Replies
15
Views
2K
  • Programming and Computer Science
Replies
2
Views
3K
  • Programming and Computer Science
Replies
14
Views
4K
  • Advanced Physics Homework Help
Replies
3
Views
2K
  • Programming and Computer Science
Replies
7
Views
6K
Back
Top