N-body problem: Difference between revisions
m
syntax highlighting fixup automation
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
|||
Line 17:
{{trans|Python}}
<
T NBody
Line 91:
print("\nCycle #.".format(i + 1))
nb.simulate()
nb.printResults()</
{{out}}
Line 154:
===The Implementations===
====Textual or Console====
<syntaxhighlight lang="c">
#include<stdlib.h>
#include<stdio.h>
Line 276:
return 0;
}
</syntaxhighlight>
Input file: same data as the TCL example.
<pre>
Line 414:
</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.
<syntaxhighlight lang="c">
#include<graphics.h>
#include<stdlib.h>
Line 542:
return 0;
}
</syntaxhighlight>
===Further work===
Line 549:
=={{header|C sharp|C#}}==
{{trans|D}}
<
using System.IO;
Line 698:
}
}
}</
{{out}}
<pre>Contents of nbody.txt
Line 816:
=={{header|C++}}==
{{trans|Java}}
<
#include <iomanip>
#include <iostream>
Line 990:
return 0;
}</
{{out}}
<pre>Contents of nbody.txt
Line 1,108:
=={{header|D}}==
{{trans|Kotlin}}
<
import std.array;
import std.conv;
Line 1,284:
nb.printResults();
}
}</
{{out}}
<pre>Contents of nbody.txt
Line 1,452:
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...
<syntaxhighlight lang="fortran">
SUBROUTINE PLOTS(DX,DY,DEVICE)
INTEGER DEVICE(10)
Line 2,071:
9001 FORMAT (' Closest approach:',E15.6)
END
</syntaxhighlight>
===Results===
Line 2,190:
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.
<syntaxhighlight lang="fortran">
Calculate some parameters for the solar system of Sun, Earth and Moon.
IMPLICIT REAL*8 (A-Z) !No integers need apply.
Line 2,260:
WRITE (6,*) 2*PI,"2Pi: distance per year for R = 1."
END
</syntaxhighlight>
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>
Line 2,435:
=={{header|Go}}==
{{trans|C}}
<
import (
Line 2,556:
}
}
}</
Contents of nbody.txt:
Line 2,697:
<li>Velocity is { vx, vy, vz }.</li>
</ul>
<
I =: 0&{"1
M =: 1&{"1
Line 2,730:
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'
Line 2,742:
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
Line 2,756:
(i.y),.(y#m),.p,.v
)
</syntaxhighlight>
<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:
Line 2,781:
¦ ¦2 1 _0.5 _0.866025 0 0.658037 _0.379918 0 ¦
+------------------------------------------------------------------------------+
</syntaxhighlight>
<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>
<syntaxhighlight lang="j">
v =: (1%3)^(1%4)
dt =: 0.001
Line 2,801:
¦ 0.080507 _0.0400152 0¦
+------------------------+
</syntaxhighlight>
<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
Line 2,864:
+--+--------+-------+---------¦
¦19¦3.04673 ¦2.06227¦0.0916438¦
+-----------------------------+</
<h5>Sensitivity to discretization</h5>
<p>Although the static-equilibrium configuration is stable, its
Line 2,870:
can examine the sensitivity to time increment, taking position error after
one rotation as a stability metric:</p>
<
dt =: y
GENDT0 3
Line 2,896:
+------+----------¦
¦1 ¦6.0354 ¦
+-----------------+</
=={{header|Java}}==
{{trans|Kotlin}}
<
import java.nio.file.Files;
import java.nio.file.Path;
Line 3,049:
}
}
}</
{{out}}
<pre>Contents of nbody.txt
Line 3,167:
=={{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),
Line 3,188:
simresult = nbodysim(bodies, timespan)
plot(simresult)
</syntaxhighlight>
=={{header|K}}==
Line 3,194:
<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
sc:{(sin 2*pi*x%y;-cos 2*pi*x%y)} / velocities</
<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
t::0;dt::0.01;rot:2 / time, time increment, rotational period</
<h4>dynamic equilibrium</h4>
<p>Within a narrow range of initial velocities the system oscillates between more and less kinetic and potential energy.</p>
<
/ v*::0.7;rot::2</
<h4>unstable</h4>
<p>The slightest perturbation from symmetry causes the system to become unstable.</p>
<
/ v+::(0 0;0 0;0 0.0001);rot:4 / wait for it ... wait for it ...</
<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
Line 3,218:
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
Line 3,230:
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
Line 3,242:
cs:{(cos 2*pi*x%y;sin 2*pi*x%y;0)} / z-axis
sc:{(sin 2*pi*x%y;-cos 2*pi*x%y;0)} / objects on screen go round and round</
=={{header|Kotlin}}==
{{trans|C}}
<
import java.io.File
Line 3,367:
nb.printResults()
}
}</
{{out}}
Line 3,487:
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
"InverseSquare", {<|"Mass" -> 1, "Position" -> {0, 0},
"Velocity" -> {0, .5}|>,
Line 3,493:
<|"Mass" -> 1, "Position" -> {0, 1}, "Velocity" -> {0, 0}|>}, 5];
{"Time:", N[#], "Positions: ", data[All, "Position", #], "Velocities: ", data[All, "Velocity", #]} & /@ Subdivide[0, 5, 20] // Grid
ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 4}]</
{{out}}
<pre>Time: 0. Positions: {{0.,0.},{1.,1.},{0.,1.}} Velocities: {{0.,0.5},{0.,-0.5},{0.,0.}}
Line 3,521:
=={{header|Nim}}==
{{trans|Go}}
<
type Vector = tuple[x, y, z: float]
Line 3,623:
echo "\nCycle ", step
sim.step()
sim.printResults()</
{{out}}
Line 3,733:
We'll show only the positions in the output.
<
global au = 150e9; # astronomical unit
Line 3,775:
disp(lsode('ABCdot',ABC, t)(1:20,1:9));
</syntaxhighlight>
{{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
Line 3,800:
=={{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.
<syntaxhighlight lang="pascal">
{$N+ Crunch only the better sort of numbers, thanks.}
{$M 50000,24000,24000] {Stack,minheap,maxheap.}
Line 4,648:
TextMode(AsItWas);
END.
</syntaxhighlight>
===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.
<syntaxhighlight lang="pascal">
{$N+ Crunch only the better sort of numbers, thanks.}
Program Swirl; Uses Graph, Crt;
Line 4,908:
WriteLn('Calculation steps:',nstep);
END.
</syntaxhighlight>
=={{header|Perl}}==
Cycle through one Neptunian year.
<
use warnings;
Line 4,974:
display($t,\@xs, \@ys) unless $t % 1000;
}
display($steps, \@xs, \@ys);</
{{out}}
<pre> Jupiter Saturn Uranus Neptune
Line 4,999:
=={{header|Phix}}==
{{trans|Kotlin}}
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">origin</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">},</span>
Line 5,083:
<span style="color: #000000;">printResults</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</
Output same as Python (and C and C# and D and Kotlin and Go)
=={{header|Python}}==
<
class Vector:
Line 5,197:
print "\nCycle %d" % (i + 1)
nb.simulate()
nb.printResults()</
{{out}}
<pre>Contents of nbody.txt
Line 5,328:
To keep things compact, we'll only display the first 20 lines of output.
<syntaxhighlight lang="raku"
multi infix:<+>(@a, @b) { @a Z+ @b }
multi infix:<->(@a, @b) { @a Z- @b }
Line 5,392:
) {
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
Line 5,417:
=={{header|Swift}}==
<
public struct Vector {
Line 5,596:
print()
}
</syntaxhighlight>
{{out}}
Line 5,719:
=={{header|Tcl}}==
{{works with|Tcl|8.6}}
<
set G 0.01
Line 5,773:
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}}
simulate $m $p $v 20</
{{out}}
<pre>
Line 5,807:
{{libheader|Wren-ioutil}}
{{libheader|Wren-fmt}}
<
import "/fmt" for Fmt
Line 5,926:
nb.simulate()
nb.printResults()
}</
{{out}}
|