Gravity Simulator
http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl
General >> Discussion >> Runge Kutta 4th order integrator
http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl?num=1306976185

Message started by Tony on 06/01/11 at 17:56:09

Title: Runge Kutta 4th order integrator
Post by Tony on 06/01/11 at 17:56:09

Try the new beta version available at http://orbitsimulator.com/gravity/beta/GravitySimulatorBeta2011June01.exe.  
There is no .dll file to download with this one, but soon I'll add one to speed it up even more.  Any comments about the new Runge Kutta routine should be made in this thread so the "Beta" thread remains more or less a log of new releases.

To use the RK4 routine, simply select it from the "Beta" menu.

Title: Re: Runge Kutta 4th order integrator
Post by wil on 06/11/11 at 17:06:46

I see You uses multiple arrays:

dx = objx(j) - objx(k)
dy = objy(j) - objy(k)
dz = objz(j) - objz(k)
...

Use one array of structures, and pointers (or references) in loops:

TObj *a = obj[i], *b = obj[j];

dx = a->x - b->x;
...

it will be faster, because offsets are calculated only ones, not eight or more.

Title: Re: Runge Kutta 4th order integrator
Post by Tony on 10/30/11 at 18:42:52

Sorry, Wil.  I didn't mean to ignore your post.  When I get back into programming mode, I'll try your suggestion.

The lastest beta version contains 4 integrators:  Euler-Cromer, Runge Kutta 4, Leap frog, and Verlet.  Here's the source code for Verlet and Leap frog:

Verlet:

Code:
   If Not firstTimeThrough Then
       For k = 1 To objects
           a1x(k) = 0
           a1y(k) = 0
           a1z(k) = 0
       Next k

       For k = 1 To objects
           Mu = ObjMass(k) * G
           For j = k + 1 To objects
               Mu2 = ObjMass(j) * G
               dx = (objx(j) - objx(k))
               dy = (objy(j) - objy(k))
               dz = (objz(j) - objz(k))
               D = Sqr(dx * dx + dy * dy + dz * dz)
               If Mu > 0 Then 'ignore massless particles
                   a1 = -(1 / D) * (1 / D) * Mu
                   a1x(j) = a1x(j) + (dx / D) * a1
                   a1y(j) = a1y(j) + (dy / D) * a1
                   a1z(j) = a1z(j) + (dz / D) * a1
               End If

               If Mu2 > 0 Then 'ignore massless particles
                   a1 = -(1 / D) * (1 / D) * Mu2
                   a1x(k) = a1x(k) - (dx / D) * a1
                   a1y(k) = a1y(k) - (dy / D) * a1
                   a1z(k) = a1z(k) - (dz / D) * a1
               End If
           Next j
       Next k
       StepEven = False
       firstTimeThrough = True
   End If

   If StepEven Then
       For k = 1 To objects
           a1x(k) = 0
           a1y(k) = 0
           a1z(k) = 0
       Next k

       For k = 1 To objects
           objx(k) = objx(k) + objvx(k) * timescale + 0.5 * a2x(k) * timescale * timescale
           objy(k) = objy(k) + objvy(k) * timescale + 0.5 * a2y(k) * timescale * timescale
           objz(k) = objz(k) + objvz(k) * timescale + 0.5 * a2z(k) * timescale * timescale
       Next k

       For k = 1 To objects
           Mu = ObjMass(k) * G
           For j = k + 1 To objects
               Mu2 = ObjMass(j) * G
               dx = (objx(j) - objx(k))
               dy = (objy(j) - objy(k))
               dz = (objz(j) - objz(k))
               D = Sqr(dx * dx + dy * dy + dz * dz)
               If Mu > 0 Then 'ignore massless particles
                   a1 = -(1 / D) * (1 / D) * Mu
                   a1x(j) = a1x(j) + (dx / D) * a1
                   a1y(j) = a1y(j) + (dy / D) * a1
                   a1z(j) = a1z(j) + (dz / D) * a1
               End If

               If Mu2 > 0 Then 'ignore massless particles
                   a1 = -(1 / D) * (1 / D) * Mu2
                   a1x(k) = a1x(k) - (dx / D) * a1
                   a1y(k) = a1y(k) - (dy / D) * a1
                   a1z(k) = a1z(k) - (dz / D) * a1
               End If
           Next j
       Next k

       For k = 1 To objects
           objvx(k) = objvx(k) + 0.5 * (a1x(k) + a2x(k)) * timescale
           objvy(k) = objvy(k) + 0.5 * (a1y(k) + a2y(k)) * timescale
           objvz(k) = objvz(k) + 0.5 * (a1z(k) + a2z(k)) * timescale
       Next k

   Else
       For k = 1 To objects
           a2x(k) = 0
           a2y(k) = 0
           a2z(k) = 0
       Next k

       For k = 1 To objects
           objx(k) = objx(k) + objvx(k) * timescale + 0.5 * a1x(k) * timescale * timescale
           objy(k) = objy(k) + objvy(k) * timescale + 0.5 * a1y(k) * timescale * timescale
           objz(k) = objz(k) + objvz(k) * timescale + 0.5 * a1z(k) * timescale * timescale
       Next k

       For k = 1 To objects
           Mu = ObjMass(k) * G
           For j = k + 1 To objects
               Mu2 = ObjMass(j) * G
               dx = (objx(j) - objx(k))
               dy = (objy(j) - objy(k))
               dz = (objz(j) - objz(k))
               D = Sqr(dx * dx + dy * dy + dz * dz)
               If Mu > 0 Then 'ignore massless particles
                   a2 = -(1 / D) * (1 / D) * Mu
                   a2x(j) = a2x(j) + (dx / D) * a2
                   a2y(j) = a2y(j) + (dy / D) * a2
                   a2z(j) = a2z(j) + (dz / D) * a2
               End If

               If Mu2 > 0 Then 'ignore massless particles
                   a2 = -(1 / D) * (1 / D) * Mu2
                   a2x(k) = a2x(k) - (dx / D) * a2
                   a2y(k) = a2y(k) - (dy / D) * a2
                   a2z(k) = a2z(k) - (dz / D) * a2
               End If
           Next j
       Next k

       For k = 1 To objects
           objvx(k) = objvx(k) + 0.5 * (a1x(k) + a2x(k)) * timescale
           objvy(k) = objvy(k) + 0.5 * (a1y(k) + a2y(k)) * timescale
           objvz(k) = objvz(k) + 0.5 * (a1z(k) + a2z(k)) * timescale
       Next k
       
   End If 'if stepeven
   StepEven = Not StepEven


Leap frog:

Code:
   If Not firstTimeThrough Then
       For k = 1 To objects
           a1x(k) = 0
           a1y(k) = 0
           a1z(k) = 0
       Next k
       For k = 1 To objects
          Mu = ObjMass(k) * G
          For j = k + 1 To objects
              Mu2 = ObjMass(j) * G
              dx = (objx(j) - objx(k))
              dy = (objy(j) - objy(k))
              dz = (objz(j) - objz(k))
              D = Sqr(dx * dx + dy * dy + dz * dz)
              If Mu > 0 Then 'ignore massless particles
                  a1 = -(1 / D) * (1 / D) * Mu
                  a1x(j) = a1x(j) + (dx / D) * a1
                  a1y(j) = a1y(j) + (dy / D) * a1
                  a1z(j) = a1z(j) + (dz / D) * a1
              End If

              If Mu2 > 0 Then 'ignore massless particles
                  a1 = -(1 / D) * (1 / D) * Mu2
                  a1x(k) = a1x(k) - (dx / D) * a1
                  a1y(k) = a1y(k) - (dy / D) * a1
                  a1z(k) = a1z(k) - (dz / D) * a1
              End If
          Next j
      Next k
      StepEven = True
      firstTimeThrough = True
   End If

   If StepEven Then
       For k = 1 To objects
           a2x(k) = 0
           a2y(k) = 0
           a2z(k) = 0
       Next k

       For k = 1 To objects
           objx(k) = objx(k) + objvx(k) * timescale + timescale * timescale * a1x(k) * 0.5
           objy(k) = objy(k) + objvy(k) * timescale + timescale * timescale * a1y(k) * 0.5
           objz(k) = objz(k) + objvz(k) * timescale + timescale * timescale * a1z(k) * 0.5
       Next k

       For k = 1 To objects
           Mu = ObjMass(k) * G
           For j = k + 1 To objects
               Mu2 = ObjMass(j) * G
               dx = (objx(j) - objx(k))
               dy = (objy(j) - objy(k))
               dz = (objz(j) - objz(k))
               D = Sqr(dx * dx + dy * dy + dz * dz)
               If Mu > 0 Then 'ignore massless particles
                   a2 = -(1 / D) * (1 / D) * Mu
                   a2x(j) = a2x(j) + (dx / D) * a2
                   a2y(j) = a2y(j) + (dy / D) * a2
                   a2z(j) = a2z(j) + (dz / D) * a2
               End If
               If Mu2 > 0 Then 'ignore massless particles
                   a2 = -(1 / D) * (1 / D) * Mu2
                   a2x(k) = a2x(k) - (dx / D) * a2
                   a2y(k) = a2y(k) - (dy / D) * a2
                   a2z(k) = a2z(k) - (dz / D) * a2
               End If
           Next j
       Next k

   Else 'stepeven
       For k = 1 To objects
           a1x(k) = 0
           a1y(k) = 0
           a1z(k) = 0
       Next k

       For k = 1 To objects
           objx(k) = objx(k) + objvx(k) * timescale + timescale * timescale * a2x(k) * 0.5
           objy(k) = objy(k) + objvy(k) * timescale + timescale * timescale * a2y(k) * 0.5
           objz(k) = objz(k) + objvz(k) * timescale + timescale * timescale * a2z(k) * 0.5
       Next k

       For k = 1 To objects
           Mu = ObjMass(k) * G
           For j = k + 1 To objects
               Mu2 = ObjMass(j) * G
               dx = (objx(j) - objx(k))
               dy = (objy(j) - objy(k))
               dz = (objz(j) - objz(k))
               D = Sqr(dx * dx + dy * dy + dz * dz)
               If Mu > 0 Then 'ignore massless particles
                   a1 = -(1 / D) * (1 / D) * Mu
                   a1x(j) = a1x(j) + (dx / D) * a1
                   a1y(j) = a1y(j) + (dy / D) * a1
                   a1z(j) = a1z(j) + (dz / D) * a1
               End If
               If Mu2 > 0 Then 'ignore massless particles
                   a1 = -(1 / D) * (1 / D) * Mu2
                   a1x(k) = a1x(k) - (dx / D) * a1
                   a1y(k) = a1y(k) - (dy / D) * a1
                   a1z(k) = a1z(k) - (dz / D) * a1
               End If
           Next j
       Next k


   End If ' stepeven
       For k = 1 To objects
           objvx(k) = objvx(k) + timescale * (a1x(k) + a2x(k)) * 0.5
           objvy(k) = objvy(k) + timescale * (a1y(k) + a2y(k)) * 0.5
           objvz(k) = objvz(k) + timescale * (a1z(k) + a2z(k)) * 0.5
       Next k

   StepEven = Not StepEven

Title: Re: Runge Kutta 4th order integrator
Post by wil on 11/08/11 at 17:47:05

Too many divisions.
div is very slow (like sqrt) - about 40 cycles.

r = dx * dx + dy * dy + dz * dz

iD3 = 1.0/(Sqrt(r)*r)       // = 1 / D^3/2

an now:

If Mu > 0 Then 'ignore massless particles
  a1 = -Mu  
  a1x(j) = a1x(j) + dx*iD3 * a1
...


Only one division instead of... 12!

Gravity Simulator » Powered by YaBB 2.1!
YaBB © 2000-2005. All Rights Reserved.