In this chapter three concepts are explained:
The second concept is used in order to calculate a position in the past. The third concept is used to find the maximum position of Mercury.
In order to simulate Newtons Law three parameters are calculated:
acceleration : a, ax and ay
speed : v, vx and vy
position : s, x and y
At each iteration or simulation cycle those three parameters are more or less calculated as follows:
an = function of m1, m2 and distance
vn = vn + an * t
sn = sn + vn * t
The problem with this method is accuracy. When the initial values are zero this is not a problem. On the other hand when the initial speed is 100 and acceleration small then after each iteration the speed stays more or less constant but the position not. The position increases linear with at each iteration becoming less accurate because the influence of a is decimated.
Two improve accuracy in x direction 3 accumulators are defined for each of the objects used: two for acceleration (Sa and SSa) and one for speed (Sv). The same for y direction
As a result for each iteration only the values of the accumulators have to be calculated and not speed and position. The position is only necessary if you want to display the place an object. The accuracy of the distance between two objects (necessary for acceleration) is drastically improved.
Initial conditions
m0 = mass of object 0
m1 = mass of object 1
s0 = distance between object 0 and object 1 at t = 0 or iteration 0
v0 = speed of object 0 at t = 0 or iteration 0
t = delta t = time period between each iteration (calculation cycle)
At iteration 1 the following values are calculated:
a0 = m1 * G / r² r² = s0 * s0
a0x = a0 * dx / s0
v1 = v0 + a0 * t
s1 = s0 + v1 * t = s0 + v0 * t + a0 * t² t² = t * t
At iteration 2 the following values are calculated:
a1 = m1 * G / r² r² = s1 * s1
a1x = a1 * dx / s1
v2 = v1 + a1 * t = v0 + (a0 + a1) * t
s2 = s1 + v2 * t = s0 + v0 * t + a0 * t² + v0 * t + (a0 + a1) * t²
= s0 + 2 * v0 * t + (2 * a0 + a1) * t²
At iteration 3 the following values are calculated:
a2 = m1 * G / r² r² = s2 * s2
a2x = a2 * G / s2
v3 = v2 + a2 * t = v0 + (a0 + a1 + a2) * t
s3 = s2 + v3 * t =
= s0 + 3 * v0 * t + (3 * a0 + 2 * a1 + a0) * t
At iteration n the following values are calculated:
an = m1 * G / r² r² = sn-1 * sn-1
anx = an * dx / sn-1
any = an * dy / sn-1
San = a0 + a1 + a2 .... an = San-1 + an
SSan = Sa0 + Sa1 + Sa2 .... San = SSan-1 + San
Svn = n * v0 = Svn-1 + v0
vnt = v0 + San * t
snt = s0 + Svn * t + SSan * t²
for x and y the following applies
vxnt = vx0 + Sanx * t (2.1)
vynt = vy0 + Sany * t (2.2)
sxnt = x0 + Svnx * t + SSanx * t² (2.3)
synt = y0 + Svny * t + SSany * t² (2.4)
Sanx = a0x + a1x + a2x .... anx = San-1x + anx (2.5)
SSanx = Sa0x + Sa1x + Sa2x .... Sanx = SSan-1x + Sanx (2.6)
Svnx = n * v0x = Svn-1x + v0x (2.7)
Sany = a0y + a1y + a2y .... any = San-1y + any (2.8)
SSany = Sa0y + Sa1y + Sa2y .... Sany = SSan-1y + Sany (2.9)
Svny = n * v0y = Svn-1y + v0y (2.10)
For each iteration only the equations 2.5 until 2.10 have to be calculated.
Using equation 2.3 the distance between two objects 1 and 2 is calculated as follows:
s = (x01 - x02) + (Svnx1 - Svnx2) * t + (SSanx1 - SSanx2) * t² (2.11)
The brackets are important. First all the values in between the brackets are calculated and then the sum is calculated thereby minimising the influence of x0
Two values are measured:
x0 at iteration = 0
x1 at iteration = 1
What is the value xn for 0
xn = a + bn
x0 = a for n = 0
x1 = a + b for n = 1
b = x1 - x0
xn = x0 + (x1 - x0) * n
for yn the following applies:
yn = y0 + (y1 - y0) * n
Three values are measured:
x0 at n = 0
x1 at n = 1
x2 at n = 2
What is the value xn for 0
xn = a + bn + cn²
x0 = a for n = 0
x1 = a + b + c for n = 1
x2 = a + 2b + 4c for n = 2
x0 = a
x1 = x0 + b + c * -2
x2 = x0 + 2b + 4c * 1
x2 - 2x1 = -x0 + 2c
c = x0/2 - x1 + x2/2
b = x1 - x0 - c = x1 - x0 - x0/2 + x1 - x2/2
= -3x0/2 + 2x1 -x2/2
xn = x0 + (-3x0 + 4x1 - x2) * n / 2 + (x0 - 2x1 + x2) * n² / 2 (2.12)
For yn the following applies:
yn = y0 + (-3y0 + 4y1 - y2) * n / 2 + (y0 - 2y1 + y2) * n² / 2 (2.13)
Example:
x0 = x0
x1 = x0 + v1.t
x2 = x1 + v2.t = x0 + v2.t + v1.t
xnt = x0 + (-3x0 + 4x0 + 4v1.t - x0 - v2.t - v1.t) * n / 2
+ ( x0 - 2x0 - 2v1.t + x0 + v2.t + v1.t) * n² / 2
xnt = x0 + (3v1.t - v2.t ) * n / 2 + (- v1.t + v2.t ) * n² / 2
Three values are measured:
xn at t = n
xn1 at t = n-1
xn2 at t = n-2
What is the value xnt for n
xnt = xn + (-3xn + 4xn1 - xn2) * n/2 + (xn - 2xn1 + xn2) * n²/2
xn = x0 + Svx.t + SSax.t²
xn-1 = x0 + Svx-1.t + SSax-1.t²
xn-2 = x0 + Svx-2.t + SSax-2.t²
xnt = x0 + Svx.t + SSax.t²
+{ -3x0 -3Svx.t -3SSax.t² +4x0 +4Svx-1.t +4SSax-1.t²
- x0 - Svx-2.t - SSax-2.t² } * n/2
+{ x0 + Svx.t + SSax.t² -2x0 -2Svx-1.t -2SSax-1.t²
+ x0 + Svx-2.t + SSax-2.t² } * n²/2
xnt = x0 + Svx.t + SSax.t²
+{ -3Svx.t -3SSax.t² +4Svx-1.t +4SSax-1.t²
- Svx-2.t - SSax-2.t² } * n/2
+{ Svx.t + SSax.t² -2Svx-1.t -2SSax-1.t²
+ Svx-2.t + SSax-2.t² } * n²/2
Svx = n * vx0
Svx-1 = (n-1) * vx0
Svx-2 = (n-2) * vx0
-3Svx + 4Svx-1 - Svx-2 = -3n*vx0 +4n*vx0 -4vx0 -n*vx0 +2vx0 = -2vx0
Svx - 2Svx-1 + Svx-2 = n*vx0 -2n*vx0 +2vx0 +n*vx0 -2vx0 = 0
xnt = x0 + Svx.t + SSax.t²
+{ -2vx0.t -3SSax.t² +4SSax-1.t² - SSax-2.t² } * n/2
+{ SSax.t² -2SSax-1.t² + SSax-2.t² } * n²/2
SSax = SSax-1 + Sax
SSax-1 = SSax-2 + Sax-1
SSax = SSax-2 + Sax-1 + Sax
xnt = x0 + Svx.t + SSax.t²
+{ -2vx0.t - 3SSax-2.t² -3Sax-1.t² -3Sax.t²
+4SSax-2.t² +4Sax-1.t² - SSax-2.t² } * n/2
+{ + SSax-2.t² + Sax-1.t² + Sax.t²
-2SSax-2.t² -2Sax-1.t² + SSax-2.t² } * n²/2
Remove SSax-2.t²
xnt = x0 + Svx.t + SSax.t²
+{ -2vx0.t -3Sax-1.t² -3Sax.t² +4Sax-1.t²} * n/2
+{ + Sax-1.t² + Sax.t² -2Sax-1.t²} * n²/2
xnt = x0 + Svx.t + SSax.t²
+{ -2vx0.t + Sax-1.t² -3Sax.t²} * n/2
+{ - Sax-1.t² + Sax.t²} * n²/2
Equation in y direction is identical.
xnt (n=0) = x0 + Svx.t + SSax.t²
xn-1t (n=1) = x0 + Svx.t + SSax.t² -vx0.t -Sax.t²
xn-2t (n=2) = x0 + Svx.t + SSax.t² -2vx0.t -Sax.t² -Sax-1.t²
For a function :
f(t) = a + b * t + c * t²
the maximum is at :
df/dt = 0
b + 2c * t = 0
t = - b / 2c
in equation 2.12 the values of of a, b and c are:
s0 = sqrt (x0²+ y0²)
etc for s1 and s2
a = s0
b = (-3s0 + 4s1 - s2)/2
c = (s0 - 2s1 + s2) / 2
Using those values the maximum is at :
3s0 - 4s1 + s2
t = ---------------
2s0 - 4s1 + 2s2
Return back to INDEX.TXT