09/22/19 at 08:01:06 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.

 Pages: 1
 Runge Kutta 4th order integrator (Read 6670 times)
 wil Uploader I Love YaBB 2! Posts: 41 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 IP Logged
 Tony YaBB Administrator Posts: 1060 Gender: 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 IP Logged
 wil Uploader I Love YaBB 2! Posts: 41 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 IP Logged
 Pages: 1
 Forum Jump: ----------------------------- General -----------------------------=> Discussion