Symmetrical n-body problem

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

Symmetrical n-body problem

Post by Lothar Schirm »

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

Post by paul doe »

Somebody has been watching too much 'The problem of the three bodies' as of lately, no? :D
(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

Post by BasicCoder2 »

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

Post by Lothar Schirm »

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
Post Reply