## Symmetrical n-body problem

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Lothar Schirm
Posts: 443
Joined: Sep 28, 2013 15:08
Location: Germany

### Symmetrical n-body problem

Long time ago, a code for simulation of the n-body system has been posted in this forum: viewtopic.php?t=10646. On Youtube, some videos show simulations of a special case of the n-body problem e.g. https://www.youtube.com/watch?v=D2YhKaANbWE or https://www.youtube.com/watch?v=Lbkr5C1i4Uo: a number of bodies with the same mass is located with equal distances on a circle, and they begin to move with the same speed along a tangent to this circle. Lagrange examined such a system of three bodies already in 1772. Though such a system is fully symmetric, it is not stable. In numeric simulations it begins to become chaotic after a while due to small deviations (a numeric simulation is never exact).

I used the Runge-Kutta method for simulation. You can try different number of bodies, initial velocities and time interval between two positions.
Code:

Code: Select all

``````'===============================================================================
' Simulation of a symmetric n-body problem
' NBody_symmetric.bas
' Created on May, 02 2024
'===============================================================================

Dim Shared As Integer n = 2, k
Dim Shared As Double v0 = 0.76, xymax = 3.0, dt = 1e-5, _
colour(), x0(), y0(), vx0(), vy0(), x1(), y1(), vx1(), vy1()
Dim As String key

Function Edit(ByVal text As String) As String
'Simple line editor to edit a predefined text (only BACKSPACE And ENTER are supported)

Dim As Integer row, col
Dim As String s, key

row = CsrLin
col = Pos
s = text

Do
Locate row, col
Print Space(Len(s) + 1)
Locate row, col
Print s;
Do
Sleep 10
key = InKey
Loop Until Len(key) = 1
Select Case ASC(key)
Case 13
Exit Do
Case 8
s = Left(s, Len(s) - 1)
Case Else
If Asc(key) > 30 Then s = s & key
End Select
Loop

Print 'new line

Return s

End Function

Sub Console()
'Console window with information and data input

Screen 0
Width 80, 25
Shell "Title " & "Simulation of a symmetric n-body problem"

Print "Initial conditions: a number of bodies are located at equal distances on "
Print "a circle and move along a tangent to this circle with the same initial "
Print "velocity. In the simulation, the radius of the circle, all masses and the "
Print "gravity constant are set to 1. The user can select the number of bodies "
Print "(2 to 7), the initial velocity, the time interval between two calculated "
Print "positions of each body and the mathematical size of the graphic window. "
Print "The program starts with a set of parameters for 3 bodies. The initial "
Print "velocity should be increased when the number of bodies is increased."
Print

Print "Edit parameters here. You can use the BACKSPACE and the ENTER key."
Print "Number of bodies (2 to 7): ";
n = Val(Edit(Str(n + 1))) - 1
Print "Initial velocity: ";
v0 = Val(Edit(Str(v0)))
Print "Time interval: ";
dt = Val(Edit(Str(dt)))
Print "Size of graphic window +/-";
xymax = Val(Edit(Str(xymax)))

ReDim colour(n), x0(n), y0(n), vx0(n), vy0(n), x1(n), y1(n), vx1(n), vy1(n)

End Sub

Function Func(Byval i As Integer, Byval t As Double, u() As Double) As Double
'Function used by the Sub DiffSysn. For each body with index k (= 0 to n), the
'following system of differential equations is solved:
'dx / dt = vx
'dy / dt = vy
'dvx / dt = -Sum( m(j) * (x(k) - x(j)) / r(k, j)^3 ), j = 0 to n, j <> k
'dvy / dt = -Sum( m(j) * (y(k) - y(j)) / r(k, j)^3 ), j = 0 to n, j <> k
'r(j, k) = sqr((x - x0(j))^2 + (y - y0(j))^2

Dim As Double x, y, vx, vy, r, sum
Dim As Integer j

x = u(1)
y = u(2)
vx = u(3)
vy = u(4)

Select Case i
Case 1
Return vx
Case 2
Return vy
Case 3
sum = 0
For j = 0 To n
If j <> k Then
r = Sqr((x - x0(j))^2 + (y - y0(j))^2)
sum = sum - (x - x0(j)) / r^3
End If
Next
Return sum
Case 4
sum = 0
For j = 0 To n
If j <> k Then
r = Sqr((x - x0(j))^2 + (y - y0(j))^2)
sum = sum - (y - y0(j)) / r^3
End If
Next
Return sum
End Select

End Function

Sub Diffsysn(ByVal x0 As Double, y0() As Double, ByVal dx As Double, Byref x1 As Double, y1() As Double)
'Stepwise solution of a system of first-order differential equations
'dy(i)/dx = Func(i, x, y(1), ... y(n)), i = 1 to n
'(Runge-Kutta method)
'Parameters:
'- x0, y0(1) to y0(n) = initial values
'- dx = stepwith
'- x1 = x0 + dx
'- y1(1) to y1(n) = results at x1
'Required function:
'Function Func(i As Integer, x As Double, y() As Double)
'Select Case i
'Case 1: Func = ... (Functions of x, y(1) to y(n))
'Case 2: Func = ...
'....
'End Select
'End Function

Dim As Integer n, i
n = Ubound(y0)
Dim As Double k1(1 To n), k2(1 To n), k3(1 To n), k4(1 To n), z(1 To n)

x1 = x0 + dx
For i = 1 To n
k1(i) = dx * Func(i, x0, y0())
Next i

For i = 1 To n
z(i) = y0(i) + k1(i) / 2
Next i
For i = 1 To n
k2(i) = dx * Func(i, x0 + dx / 2, z())
Next i

For i = 1 To n
z(i) = y0(i) + k2(i) / 2
Next i
For i = 1 To n
k3(i) = dx * Func(i, x0 + dx / 2, z())
Next i

For i = 1 To n
z(i) = y0(i) + k3(i)
Next i
For i = 1 To n
k4(i) = dx * Func(i, x1, z())
Next i

For i = 1 To n
y1(i) = y0(i) + (k1(i) + 2 * k2(i) + 2 * k3(i) + k4(i)) / 6
Next i

End Sub

Sub Simulation()

Dim As Integer i, xp(n), yp(n)
Dim As Any Ptr img(n)
Dim As Double pi = 4 * Atn(1), phi, u0(4), u1(4), t0 = 0, t1

'Open graphic window :
Screenres 600,650
Windowtitle "Simulation"
Width 600\8, 650\16
Locate Hiword(Width), 10
Print "Stop: Press any key";

'Image buffers for bodies:
For i = 0 To n
colour(i) = 9 + i
img(i) = Imagecreate(20, 20, 0)
Circle img(i), (10, 10), 9, colour(i),,,,F
Next

'Define initial conditions:
For i = 0 To n
phi = pi / 2 - 2 * i * pi / (n + 1)
x0(i) = Cos(phi)
y0(i) = Sin(phi)
vx0(i) = v0 * Cos(phi - pi / 2)
vy0(i) = v0 * Sin(phi - pi / 2)
Next

View (0, 0) - (599, 599),, 3
Window (- xymax, - xymax) - ( xymax,  xymax)

Do

'Draw bodies:
For k = 0 To n
xp(k) = Pmap(x0(k), 0) - 10
yp(k) = Pmap(y0(k), 1) - 5
Put (Pmap(xp(k), 2), Pmap(yp(k), 3)), img(k), XOR
Next

'Calculate next positions:
For k = 0 To n
u0(1) = x0(k)
u0(2) = y0(k)
u0(3) = vx0(k)
u0(4) = vy0(k)
DiffSysn(t0, u0(), dt, t1, u1())
x1(k) = u1(1)
y1(k) = u1(2)
vx1(k) = u1(3)
vy1(k) = u1(4)
Next

'Delete bodies and draw actual position as a point:
For k = 0 To n
Put (Pmap(xp(k), 2), Pmap(yp(k), 3)), img(k), Xor
PSet (x0(k), y0(k)), colour(k)
Next

'The new positions are used as the next starting points:
t0 = t1
For k = 0 To n
x0(k) = x1(k)
y0(k) = y1(k)
vx0(k) = vx1(k)
vy0(k) = vy1(k)
Next

Loop Until Inkey <> ""

For k = 0 To n
ImageDestroy img(k)
Next

End Sub

'Main:

Do
Console()
Simulation()
Locate Hiword(Width), 10
Print Space(50);
Locate HiWord(Width), 10
Input; "Again y/n: ", key
Loop Until LCase(key) = "n"

End
``````
paul doe
Moderator
Posts: 1745
Joined: Jul 25, 2017 17:22
Location: Argentina

### Re: Symmetrical n-body problem

Somebody has been watching too much 'The problem of the three bodies' as of lately, no?
(Not that I blame you, it's amazing)

This is a fascinating problem, indeed. Very nice.
BasicCoder2
Posts: 3917
Joined: Jan 01, 2009 7:03
Location: Australia

### Re: Symmetrical n-body problem

Interesting simulation although I don't understand the math.
In numeric simulations it begins to become chaotic after a while due to small deviations (a numeric simulation is never exact).
Isn't the real world the same? We can't measure or represent a real world state variable to an infinitely exact value and thus the simulation to predict the future world states will eventually diverge from the predicted next state. Chaos theory?
Lothar Schirm
Posts: 443
Joined: Sep 28, 2013 15:08
Location: Germany

### Re: Symmetrical n-body problem

Thank you, paul doe, and best regards!

Hi BasicCoder2! Yes, the math is a little bit complex. If you search for informations on the Runge-Kutta method, you will mostly find explanations how to solve first-oder differential equations, not systems. For systems of first-order differential equations, I found a code for QuickBASIC many years ago which I modified and use now for FreeBasic. Regarding the gravitational forces between the different bodies in a two-dimensional problem, you must use the vector components in the directions of x and y, etc., so you need some knowledge about vectors and Newton laws of motion ... I studied physics and I am interested in astronomy, so writing code in this way is very interesting for me. I posted some of my works on FreeBasic Portal, using my "Simple GUI" as user interface, but in german, sorry! (https://www.freebasic-portal.de/downloa ... n-361.html.

Regarding your consideration about the possibility to predict the real world, I agree. Even in physics, which seems to be the most exact natural science, we always work with models which are simplified so far that they are able to predict something, but we have always to check that the results are in accordance with the real world.

Best regards to the other side of the world!
Lothar