N-body problem: Difference between revisions

;Task

Write a simulation of three masses interacting under gravitation and show their evolution over a number of time-steps (at least 20).

 

=={{header|C}}==

===The Implementations===

====Textual or Console====

#include<stdlib.h>

#include<stdio.h>

return 0;

}

Input file: same data as the TCL example.

<pre>

</pre>

The time step although present in the file, and read, is not used here, the loop goes on indefinitely. Requires the [http://www.cs.colorado.edu/~main/bgi/cs1300/ WinBGIm] library.

#include<graphics.h>

#include<stdlib.h>

return 0;

}

 

===Further work===

=={{header|C sharp|C#}}==

{{trans|D}}

using System.IO;

 

}

}

{{out}}

<pre>Contents of nbody.txt

=={{header|C++}}==

{{trans|Java}}

#include <iomanip>

#include <iostream>

 

return 0;

{{out}}

<pre>Contents of nbody.txt

=={{header|D}}==

{{trans|Kotlin}}

import std.array;

import std.conv;

nb.printResults();

}

{{out}}

<pre>Contents of nbody.txt

 

More seriously, later Fortran has additional routines, and confusion can arise: thus subroutine FREE was renamed to SFREE, and array INDEX renamed to INDEXS, though calling it SINDEX was the first idea - but that would require additional declarations so as to keep it integer. The habit of not bothering with minor declarations lingers on with exploratory programmes...

SUBROUTINE PLOTS(DX,DY,DEVICE)

INTEGER DEVICE(10)

9001 FORMAT (' Closest approach:',E15.6)

END

 

===Results===

 

The idea was to start with an initial state having the Sun, Earth and Moon all in a straight line along the x-axis, but after some ad-hoc messing about, confusions escalated to the degree that an ad-hoc calculation prog. was in order. This turns out to need a cube root function, and a proper solution for this raises again the utility of palindromic functions. Later Fortran supplies intrinsic functions such as EXPONENT(x) which returns the exponent part of the floating-point number, ''x'' and extracting this would be useful in devising an initial value for an iterative refinement calculation. Clearly, one wants power/3 for this, and so being able to write something like <code>EXPONENT(t) = EXPONENT(x)/3</code> would help, along with similar usage of the FRACTION(x) function - omitting details such as the remainder when the power is divided by three. The same ideas arise with the square root function, though here, SQRT is already supplied. Alas, only the SUBSTR intrinsic function of pl/i has palindromic usage.

Calculate some parameters for the solar system of Sun, Earth and Moon.

IMPLICIT REAL*8 (A-Z) !No integers need apply.

WRITE (6,*) 2*PI,"2Pi: distance per year for R = 1."

END

Notably, wanting a circular orbit to ease inspection of the results meant that the actual semi-major axis of the earth's elliptical orbit could not be used as in the first trials. Instead, the orbital period was taken as a year (wrongly so, as this includes the precession of the earth's axis of revolution, with a period of about 26,000 years) and the radius of a circular orbit having that period determined. This however is affected by whether the calculation is for the sun alone (so, a massless Earth), or, the sun and the earth together, or the sun, earth and moon together. As well, there is no z-action: the moon's orbital plane is deemed the same as that of the earth and it isn't. Results:

<pre>

=={{header|Go}}==

{{trans|C}}

 

import (

}

}

 

Contents of nbody.txt:

<li>Velocity is { vx, vy, vz }.</li>

</ul>

I =: 0&{"1

M =: 1&{"1

out =: out,D z

z

<h4>Integrator</h4>

<p>Iterate over a time interval. Plot the results (x and y coordinates over time). </p>

maxn =: 20000

require'plot'

dt NEXT^:x y

plot 0 1 { |: out

<h4>Configurations</h4>

<h5>Configuration generator</h5>

<p>Generate a radially symmetrical configuration for a given number of bodies and an initial rotational velocity</p>

GEN =: 3 : 0

1 GEN y

(i.y),.(y#m),.p,.v

)

<h5>Generate a few cases.</h5>

<p>The 3-body static equilibrium case is proposed as the basis for validation, and for comparison of alternative methods.</p>

case =. 3 : 'y;".y'

r =. 0 2 $ a:

¦ ¦2 1 _0.5 _0.866025 0 0.658037 _0.379918 0 ¦

+------------------------------------------------------------------------------+

<h5>Static equilibrium case. </h5>

<p> Bodies orbit with constant kinetic energy.</p>

[<em>sinusoidal curves of constant amplitude</em>[https://commons.wikimedia.org/wiki/File:Nbody_j_2.jpg]]

<p>Determine orbital period, and compare configuration after two rotations vs. initial position. </p>

v =: (1%3)^(1%4)

dt =: 0.001

¦ 0.080507 _0.0400152 0¦

+------------------------+

<h5>Dynamic equilibrium case. </h5>

<p> Within a narrow range of initial velocities, the bodies orbit, maintaining symetry and swapping kinetic and potential energy back and forth. The ideal system is stable, but an integrating solver will not be, due to numerical precision. </p>

[<em>sinusoidal curves of varying amplitudes</em> [https://commons.wikimedia.org/wiki/File:Nbody_j_1.jpg]]

<h5>Perturbed case. </h5>

<p>Any perturbation from symmetry, however small, will eventually cause the bodies to escape. Note that the effect demonstrated here is due to the physics of the problem, independently of the numerical precision and increment size of the solver. </p>

[<em>curves start out nice then get jumbled and diverge</em> [https://commons.wikimedia.org/wiki/File:Nbody_j_3.jpg]]

<h5>n>3</h5>

<p>We can generate and execute tests for radially symmetrical configurations of any number of bodies, run the simulations, and compare the results with an explicitly determined result:</p>

f0 =. 4 : '1%((x r y)^2)' NB. force due to nth body

f =. 4 : '(x f0 y)*(1-cos 2*pi*y%x)%x r y'"0 NB. centripetal component of force

+--+--------+-------+---------¦

¦19¦3.04673 ¦2.06227¦0.0916438¦

<h5>Sensitivity to discretization</h5>

<p>Although the static-equilibrium configuration is stable, its

can examine the sensitivity to time increment, taking position error after

one rotation as a stability metric:</p>

dt =: y

GENDT0 3

+------+----------¦

¦1 ¦6.0354 ¦

 

=={{header|Java}}==

{{trans|Kotlin}}

import java.nio.file.Files;

import java.nio.file.Path;

}

}

{{out}}

<pre>Contents of nbody.txt

Body 2 : 0.426346 -0.111425 -0.681150 | 0.234987 0.006241 -0.440245

Body 3 : -1.006089 -4.103186 -3.122591 | -0.359686 -1.177995 -0.875924</pre>

 

=={{header|Julia}}==

Uses the NBodySimulator module.

 

const stablebodies = [MassBody(SVector(0.0, 1.0, 0.0), SVector( 5.775e-6, 0.0, 0.0), 2.0),

simresult = nbodysim(bodies, timespan)

plot(simresult)

 

=={{header|K}}==

<h4>configuration generator</h4>

<p>The configurations shown here are based on those in the J task.</p>

cs:{(cos 2*pi*x%y;sin 2*pi*x%y)} / positions

<h4>static equilibrium configuration</h4>

<p>The system is stable for an initial velocity for which centrifugal and gravitational forces are in balance</p>

g::1;vs::pow[1%3;1%4];v*::vs / gravitational constant, initial velocity

<h4>dynamic equilibrium</h4>

<p>Within a narrow range of initial velocities the system oscillates between more and less kinetic and potential energy.</p>

<h4>unstable</h4>

<p>The slightest perturbation from symmetry causes the system to become unstable.</p>

<h4>physics</h4>

A1:{$[x=y;0;g*m[y]*DV[x;y]%pow[DS[x;y];3]]} / nth body centripetal acceleration component

A:{+/{A1/:[x;!n]}'!n} / total centripetal acceleration

v0:v;v+::dt*A[];p+::dt*0.5*v+v0 / update velocities and positions

t+::dt;$[t>tmax;do[msg::,"[";tick::{}];0] / update time, detect end of time

<h4>events and graphics</h4>

sc:1.25;P::{w*0.5*1+(1%sc)*x} / display units

wh::(w;h);C::wh%2 / viewbox

r:(,:'P[p]-s%2),'`RGB,',:'(2#'s)#'!n / plot the bodies

r,:,(P[1 0];BW;~,/'+text@`i$msg) / plot the initial position marker

<h4>3D</h4>

<p>Although the animation is limited to two dimensions,

we can verify that the implementation really does support three,

by changing the axis of rotation in the initialization function:</p>

sc:{(0;sin 2*pi*x%y;-cos 2*pi*x%y)} / objects on screen move up and down

 

 

cs:{(cos 2*pi*x%y;sin 2*pi*x%y;0)} / z-axis

 

=={{header|Kotlin}}==

{{trans|C}}

 

import java.io.File

nb.printResults()

}

 

{{out}}

Body 3 : -1.006089 -4.103186 -3.122591 | -0.359686 -1.177995 -0.875924

</pre>

 

=={{header|Nim}}==

{{trans|Go}}

 

type Vector = tuple[x, y, z: float]

echo "\nCycle ", step

sim.step()

 

{{out}}

We'll show only the positions in the output.

 

global au = 150e9; # astronomical unit

 

 

disp(lsode('ABCdot',ABC, t)(1:20,1:9));

{{out}}

<pre> 0.0000e+00 0.0000e+00 0.0000e+00 1.5000e+11 0.0000e+00 0.0000e+00 1.4970e+11 0.0000e+00 0.0000e+00

=={{header|Pascal}}==

This was written in Turbo Pascal, prompted by an article in an astronomy magazine about the orbit of an asteroid about a planet about a sun. It expects to receive command-line parameters for the run, and wants to use the screen graphics to animate the result. Alas, the turbo .bgi protocol no longer works on modern computers that have screens with many more dots. It contains options to select different methods for computation. An interesting point is that if the central mass is fixed in position (rather than wobbling about the centre of mass, as it does), then the Trojan orbit positions are not stable.

{$N+ Crunch only the better sort of numbers, thanks.}

{$M 50000,24000,24000] {Stack,minheap,maxheap.}

TextMode(AsItWas);

END.

 

===Relativistic orbits===

Another magazine article described a calculation for orbits twisted by the gravity around black holes. Only a two-body problem, and only just, because the orbiting body's mass is ignored.

{$N+ Crunch only the better sort of numbers, thanks.}

Program Swirl; Uses Graph, Crt;

WriteLn('Calculation steps:',nstep);

END.

 

=={{header|Perl}}==

Cycle through one Neptunian year.

use warnings;

 

display($t,\@xs, \@ys) unless $t % 1000;

}

{{out}}

<pre> Jupiter Saturn Uranus Neptune

=={{header|Phix}}==

{{trans|Kotlin}}

Output same as Python (and C and C# and D and Kotlin and Go)

 

=={{header|Python}}==

 

class Vector:

print "\nCycle %d" % (i + 1)

nb.simulate()

{{out}}

<pre>Contents of nbody.txt

(formerly Perl 6)

 

 

We use a 18-dimension vector <math>ABC</math>. The first nine dimensions are the positions of the three bodies. The other nine are the velocities. This allows us to write the dynamics as a first-temporal derivative equation, since

To keep things compact, we'll only display the first 20 lines of output.

 

multi infix:<+>(@a, @b) { @a Z+ @b }

multi infix:<->(@a, @b) { @a Z- @b }

) {

printf "t = %.02f : %s\n", $t, @ABC.fmt("%+.3e");

{{out}}

<pre>t = 0.00 : +0.000e+00 +0.000e+00 +0.000e+00 +1.500e+11 +0.000e+00 +0.000e+00 +1.497e+11 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +2.987e+04 +0.000e+00 +0.000e+00 +3.090e+04 +0.000e+00

=={{header|Swift}}==

 

 

public struct Vector {

print()

}

 

{{out}}

=={{header|Tcl}}==

{{works with|Tcl|8.6}}

 

set G 0.01

puts [lmap pos $p {format (%.5f,%.5f,%.5f) {*}$pos}]

}

Demonstrating for 20 steps:

set p {{0 0 0} {1 1 0} {0 1 1}}

set v {{0.01 0 0} {0 0 0.02} {0.01 -0.01 -0.01}}

{{out}}

<pre>

{{libheader|Wren-ioutil}}

{{libheader|Wren-fmt}}

 

class Vector3D {

nb.simulate()

nb.printResults()

 

{{out}}