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 * dziD3 = 1.0/(Sqrt(r)*r)       // = 1 / D^3/2an now: If Mu > 0 Then 'ignore massless particles   a1 = -Mu     a1x(j) = a1x(j) + dx*iD3 * a1...Only one division instead of... 12!