Welcome, Guest. Please Login.
Gravity Simulator
11/22/17 at 17:57:55
News: Registration for new users has been disabled to discourage spam. If you would like to join the forum please send me an email with your desired screen name to tony at gravitysimulator dot com.
Home Help Search Login


Pages: 1
Send Topic Print
Runge Kutta 4th order integrator (Read 4815 times)
Tony
YaBB Administrator
*****




Posts: 1051
Gender: male
Runge Kutta 4th order integrator
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.
Back to top
 
 
Email View Profile WWW   IP Logged
wil
Uploader



I Love YaBB 2!

Posts: 39
Re: Runge Kutta 4th order integrator
Reply #1 - 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.
Back to top
 
 
Email View Profile   IP Logged
Tony
YaBB Administrator
*****




Posts: 1051
Gender: male
Re: Runge Kutta 4th order integrator
Reply #2 - 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
 

Back to top
 
 
Email View Profile WWW   IP Logged
wil
Uploader



I Love YaBB 2!

Posts: 39
Re: Runge Kutta 4th order integrator
Reply #3 - 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!
Back to top
 
 
Email View Profile   IP Logged
Pages: 1
Send Topic Print