- #1

Jenab

Science Advisor

- 114

- 0

## Main Question or Discussion Point

To help me in understanding this procedure, I read the fifth chapter of THE DETERMINATION OF ORBITS by A. D. Dubyago, as translated from Russian into English by the RAND Corporation. Much of the nomenclature used in this exposition is patterned after Dubyago.

For i=1 to 3

The time of observation: t(i)

The position of the sun: Xs(i), Ys(i), Zs(i)

The unit vector in the direction of the planet: a(i), b(i), c(i)

The position vectors are in local celestial coordinates.

k = 0.01720209895 (mean motion of Earth, radians per day)

A = 0.00577551833 (days required for light to travel 1 AU)

tau(1) = k {t(3)-t(2)}

tau(2) = k {t(3)-t(1)}

tau(3) = k {t(2)-t(1)}

nc(1) = tau(1)/tau(2)

nc(3) = tau(3)/tau(2)

nu(1) = tau(1) tau(3) {1 + nc(1)} / 6

nu(3) = tau(1) tau(3) {1 + nc(3)} / 6

D = a(2) { b(3) c(1) - b(1) c(3) } + b(2) { c(3) a(1) - c(1) a(3) } + c(2) { a(3) b(1) - a(1) b(3) }

For i=1 to 3

d(i) = Xs(i) { c(1) b(3) - c(3) b(1) } + Ys(i) { a(1) c(3) - a(3) c(1) } + Zs(i) { b(1) a(3) - b(3) a(1) }

For i=1 to 3

Re(i) = { Xs(i)^2 + Ys(i)^2 + Zs(i)^2 }^0.5

For i=1 to 3

q(i) = -2 { a(i) Xs(i) + b(i) Ys(i) + c(i) Zs(i) }

K0 = { d(2) - d(1) nc(1) - d(3) nc(3) } / D

L0 = { d(1) nu(1) + d(3) nu(3) } / D

Let r(2) be the distance from the sun to the planet at time 2

Let p(2) be the distance from Earth to the planet at time 2.

We must solve these two equations simultaneously for r(2) and p(2):

p(2) = K0 - L0 / r(2)^3

r(2) = { Re(2)^2 + q(2) p(2) + p(2)^2 }^0.5

This leads to an 8th degree polynomial in r(2):

f{r(2)} = r(2)^8 - { Re(2)^2 + q(2) K0 + K0^2 } r(2)^6 + L0 {q(2) + 2 K0} r(2)^3 - L0^2

f{r(2)} = 0

We will use Newton's method to get the answer.

df/dr(2) = 8 r(2)^7 - 6 { Re(2)^2 + q(2) K0 + K0^2 } r(2)^5 + 3 L0 {q(2) + 2 K0} r(2)^2

As an initial guess at r(2)...

r2(j=0) = 1.0

Then...

Repeat over j

r2(j+1) = r2(j) - f{r2(j)} / df/dr2(j)

Until abs{ r2(j+1) - r2(j) } < 1.0E-12

Assign to r(2) the converged value from the loop:

r(2) = r2(j+1)

p(2) = K0 - L0 / r(2)^3

n(1) = nc(1) + nu(1) / r(2)^3

n(3) = nc(3) + nu(3) / r(2)^3

Q1 = n(1) Xs(1) - Xs(2) + n(3) Xs(3) + a(2) p(2)

Q2 = n(1) Ys(1) - Ys(2) + n(3) Ys(3) + b(2) p(2)

Q3 = n(1) Zs(1) - Zs(2) + n(3) Zs(3) + c(2) p(2)

Q4 = a(1) - a(3) c(1) / c(3)

Q5 = b(1) - b(3) a(1) / a(3)

Q6 = b(3) - b(1) a(3) / a(1)

Q7 = a(3) - a(1) c(3) / c(1)

p(1) = { [ Q1 - a(3) Q3 / c(3) ] / [ n(1) Q4 ] + [ Q2 - b(3) Q1 / a(3) ] / [n(1) Q5] } / 2

p(3) = { [ Q2 - b(1) Q1 / a(1) ] / [ n(3) Q6 ] + [ Q1 - a(1) Q3 / c(1) ] / [n(3) Q7] } / 2

r(1) = { Re(1)^2 + q(1) p(1) + p(1)^2 }^0.5

r(3) = { Re(3)^2 + q(3) p(3) + p(3)^2 }^0.5

Light doesn't travel instantaneously. We must correct the initial times for the amount of time it took light to travel from the planet to our eyes. This is done only once, after the first approximation, but before the 2nd approximation. It will not be done between any other successive approximations.

For i=1 to 3

tc(i) = t(i) - A p(i)

Where A is the time (in days) required for light to travel 1 astronomical unit.

tau(1) = k {tc(3)-tc(2)}

tau(2) = k {tc(3)-tc(1)}

tau(3) = k {tc(2)-tc(1)}

nc(1) = tau(1)/tau(2)

nc(3) = tau(3)/tau(2)

For i=1 to 3

x(i) = a(i) p(i) - Xs(i)

y(i) = b(i) p(i) - Ys(i)

z(i) = c(i) p(i) - Zs(i)

K(1) = { 2 [ r(2) r(3) + x(2) x(3) + y(2) y(3) + z(2) z(3) ] }^0.5

K(2) = { 2 [ r(1) r(3) + x(1) x(3) + y(1) y(3) + z(1) z(3) ] }^0.5

K(3) = { 2 [ r(1) r(2) + x(1) x(2) + y(1) y(2) + z(1) z(2) ] }^0.5

h(1) = tau(1)^2 / { K(1)^2 [ K(1)/3 + ( r(2)+r(3) )/2 ] }

h(2) = tau(2)^2 / { K(2)^2 [ K(2)/3 + ( r(1)+r(3) )/2 ] }

h(3) = tau(3)^2 / { K(3)^2 [ K(3)/3 + ( r(1)+r(2) )/2 ] }

For i=1 to 3

Begin

Queasy1 = (11/9) h(i)

Queasy2 = Queasy1

Repeat

Queasy3 = Queasy2

Queasy2 = Queasy1 / ( 1 + Queasy2)

Until abs( Queasy2 - Queasy3 ) / Queasy3 < 1E-12;

Yippy(i) = 1 + (10/11) Queasy2

End

n(1) = nc(1) yippy(2) / yippy(1)

n(3) = nc(3) yippy(2) / yippy(3)

nu(1) = nc(1) r(2)^3 { yippy(2) / yippy(1) - 1 }

nu(3) = nc(3) r(2)^3 { yippy(2) / yippy(3) - 1 }

K0 = { d(2) - d(1) nc(1) - d(3) nc(3) } / D

L0 = { d(1) nu(1) + d(3) nu(3) } / D

p(2) = K0 - L0 / r(2)^3

r(2) = { Re(2)^2 + q(2) p(2) + p(2)^2 }^0.5

f{r(2)} = r(2)^8 - { Re(2)^2 + q(2) K0 + K0^2 } r(2)^6 + L0 {q(2) + 2 K0} r(2)^3 - L0^2

f{r(2)} = 0

df/dr(2) = 8 r(2)^7 - 6 { Re(2)^2 + q(2) K0 + K0^2 } r(2)^5 + 3 L0 {q(2) + 2 K0} r(2)^2

The customary initial guess at r(2)...

r2(j=0) = 1.0

Repeat over j

r2(j+1) = r2(j) - f{r2(j)} / df/dr2(j)

Until abs{ r2(j+1) - r2(j) } < 1.0E-12

Assign to r(2) the converged value from the loop:

r(2) = r2(j+1)

p(2) = K0 - L0 / r(2)^3

Q1 = n(1) Xs(1) - Xs(2) + n(3) Xs(3) + a(2) p(2)

Q2 = n(1) Ys(1) - Ys(2) + n(3) Ys(3) + b(2) p(2)

Q3 = n(1) Zs(1) - Zs(2) + n(3) Zs(3) + c(2) p(2)

Q4 = a(1) - a(3) c(1) / c(3)

Q5 = b(1) - b(3) a(1) / a(3)

Q6 = b(3) - b(1) a(3) / a(1)

Q7 = a(3) - a(1) c(3) / c(1)

p(1) = { [ Q1 - a(3) Q3 / c(3) ] / [ n(1) Q4 ] + [ Q2 - b(3) Q1 / a(3) ] / [n(1) Q5] } / 2

p(3) = { [ Q2 - b(1) Q1 / a(1) ] / [ n(3) Q6 ] + [ Q1 - a(1) Q3 / c(1) ] / [n(3) Q7] } / 2

r(1) = { Re(1)^2 + q(1) p(1) + p(1)^2 }^0.5

r(3) = { Re(3)^2 + q(3) p(3) + p(3)^2 }^0.5

Repeat this whole section (Recursive Procedure for Successive Approximations) until your values for r(i) converge.

MORE ORBIT DETERMINATION IN POSTS TO FOLLOW.

t(1) = JD 2427255.460417 = 1 July 1933, 23h 3m 0s UT

Xs(1) = -0.169709

Ys(1) = +0.919710

Zs(1) = +0.398865

RA(1) = 19.46730000 hours

DEC(1) = -13.86869444 degrees

t(2) = JD 2427283.391181 = 29 July 1933, 21h 23m 18s UT

Xs(2) = -0.600429

Ys(2) = +0.751016

Zs(2) = +0.325697

RA(2) = 19.06218056 hours

DEC(2) = -14.11902778 degrees

t(3) = JD 2427312.342083 = 27 August 1933, 20h 12m 36s UT

Xs(3) = -0.908371

Ys(3) = +0.405220

Zs(3) = +0.175716

RA(3) = 18.98696667 hours

DEC(3) = -15.24394444 degrees

a(1) = +0.363835156

b(1) = -0.900093900

c(1) = -0.239697622

a(2) = +0.266215600

b(2) = -0.932536300

c(2) = -0.243937090

a(3) = +0.246531192

b(3) = -0.932786462

c(3) = -0.262929245

tau(1) = 0.498016281

tau(2) = 0.978484047

tau(3) = 0.480467766

nc(1) = 0.508967195

nc(3) = 0.491032805

nu(1) = 0.060177805

nu(3) = 0.059462580

d(1) = 0.015443444

d(2) = 0.018648221

d(3) = 0.017700437

D = 0.001964676

Re(1) = 1.016740339

Re(2) = 1.015193850

Re(3) = 1.010058035

q(1) = 1.970356907

q(2) = 1.879285653

q(3) = 1.296252781

K0 = 1.067106795

L0 = 1.008749516

The Lagrange equation is

f(r2) = r(2)^8 - 4.174733956 r(2)^6 + 4.048615418 r(2)^3 - 1.017575585

The first derivative of the Lagrange equation is

f'(r2) = 8 r(2)^7 - 25.048403737 r(2)^5 + 12.145846255 r(2)^2

This actually takes a while to converge with Newton's method, so I'd consider going to Danby's method for faster convergence. Anyway, it converges to

r(2) = 1.898670742

p(2) = 0.919728216

n(1) = 0.517759189

n(3) = 0.499720304

Q1 = +0.303475173

Q2 = -0.930010982

Q3 = -0.255727953

Q4 = +0.139086706

Q5 = +0.476529097

Q6 = -0.322891519

Q7 = -0.152567065

p(1) = 0.884503909

p(3) = 1.110851828

r(1) = 1.886503769

r(3) = 1.922018156

tc(1) = JD 2427255.455309

tc(2) = JD 2427283.385869

tc(3) = JD 2427312.335667

tau(1) = 0.497997293

tau(2) = 0.978461559

tau(3) = 0.480464267

nc(1) = 0.508959486

nc(3) = 0.491040514

x(1) = +0.491522618

y(1) = -1.715846574

z(1) = -0.610878483

x(2) = +0.845274999

y(2) = -1.608695948

z(2) = -0.550052825

x(3) = +1.182230625

y(3) = -1.441407547

z(3) = -0.467791433

K(1) = 3.801232986

K(2) = 3.732555561

K(3) = 3.766593197

h(1) = 0.005401695

h(2) = 0.021826229

h(3) = 0.005168609

Yippy(1) = 1.005962773

Yippy(2) = 1.023636797

Yippy(3) = 1.005707071

nu(1) = 0.061204833

nu(3) = 0.059919537

n(1) = 0.517901529

n(3) = 0.499794774

K0 = 1.067097940

L0 = 1.020939404

The Lagrange equation is

f(r2) = r(2)^8 - 4.174698414 r(2)^6 + 4.097521445 r(2)^3 - 1.042317267

The first derivative of the Lagrange equation is

f'(r2) = 8 r(2)^7 - 25.048190484 r(2)^5 + 12.292564335 r(2)^2

Newton's method converges to...

r(2) = 1.896390493

p(2) = 0.917399712

[Skipping the Q's.]

p(1) = 0.882742391

p(3) = 1.107232428

r(1) = 1.884757972

r(3) = 1.918706335

You've seen how it works. You just keep taking p(i) from the end of the latest approximation and plugging it back into the top of the next approximation and repeating the math with the improved numbers. You keep repeating the recursive section until all the values for r(i) have converged.

Jerry Abbott

**Initial Data.**For i=1 to 3

The time of observation: t(i)

The position of the sun: Xs(i), Ys(i), Zs(i)

The unit vector in the direction of the planet: a(i), b(i), c(i)

The position vectors are in local celestial coordinates.

**Handy Constants.**k = 0.01720209895 (mean motion of Earth, radians per day)

A = 0.00577551833 (days required for light to travel 1 AU)

**First Approximation.**tau(1) = k {t(3)-t(2)}

tau(2) = k {t(3)-t(1)}

tau(3) = k {t(2)-t(1)}

nc(1) = tau(1)/tau(2)

nc(3) = tau(3)/tau(2)

nu(1) = tau(1) tau(3) {1 + nc(1)} / 6

nu(3) = tau(1) tau(3) {1 + nc(3)} / 6

D = a(2) { b(3) c(1) - b(1) c(3) } + b(2) { c(3) a(1) - c(1) a(3) } + c(2) { a(3) b(1) - a(1) b(3) }

For i=1 to 3

d(i) = Xs(i) { c(1) b(3) - c(3) b(1) } + Ys(i) { a(1) c(3) - a(3) c(1) } + Zs(i) { b(1) a(3) - b(3) a(1) }

For i=1 to 3

Re(i) = { Xs(i)^2 + Ys(i)^2 + Zs(i)^2 }^0.5

For i=1 to 3

q(i) = -2 { a(i) Xs(i) + b(i) Ys(i) + c(i) Zs(i) }

K0 = { d(2) - d(1) nc(1) - d(3) nc(3) } / D

L0 = { d(1) nu(1) + d(3) nu(3) } / D

Let r(2) be the distance from the sun to the planet at time 2

Let p(2) be the distance from Earth to the planet at time 2.

We must solve these two equations simultaneously for r(2) and p(2):

p(2) = K0 - L0 / r(2)^3

r(2) = { Re(2)^2 + q(2) p(2) + p(2)^2 }^0.5

This leads to an 8th degree polynomial in r(2):

f{r(2)} = r(2)^8 - { Re(2)^2 + q(2) K0 + K0^2 } r(2)^6 + L0 {q(2) + 2 K0} r(2)^3 - L0^2

f{r(2)} = 0

We will use Newton's method to get the answer.

df/dr(2) = 8 r(2)^7 - 6 { Re(2)^2 + q(2) K0 + K0^2 } r(2)^5 + 3 L0 {q(2) + 2 K0} r(2)^2

As an initial guess at r(2)...

r2(j=0) = 1.0

Then...

Repeat over j

r2(j+1) = r2(j) - f{r2(j)} / df/dr2(j)

Until abs{ r2(j+1) - r2(j) } < 1.0E-12

Assign to r(2) the converged value from the loop:

r(2) = r2(j+1)

p(2) = K0 - L0 / r(2)^3

n(1) = nc(1) + nu(1) / r(2)^3

n(3) = nc(3) + nu(3) / r(2)^3

Q1 = n(1) Xs(1) - Xs(2) + n(3) Xs(3) + a(2) p(2)

Q2 = n(1) Ys(1) - Ys(2) + n(3) Ys(3) + b(2) p(2)

Q3 = n(1) Zs(1) - Zs(2) + n(3) Zs(3) + c(2) p(2)

Q4 = a(1) - a(3) c(1) / c(3)

Q5 = b(1) - b(3) a(1) / a(3)

Q6 = b(3) - b(1) a(3) / a(1)

Q7 = a(3) - a(1) c(3) / c(1)

p(1) = { [ Q1 - a(3) Q3 / c(3) ] / [ n(1) Q4 ] + [ Q2 - b(3) Q1 / a(3) ] / [n(1) Q5] } / 2

p(3) = { [ Q2 - b(1) Q1 / a(1) ] / [ n(3) Q6 ] + [ Q1 - a(1) Q3 / c(1) ] / [n(3) Q7] } / 2

r(1) = { Re(1)^2 + q(1) p(1) + p(1)^2 }^0.5

r(3) = { Re(3)^2 + q(3) p(3) + p(3)^2 }^0.5

**Correction for Planetary Abberation.**Light doesn't travel instantaneously. We must correct the initial times for the amount of time it took light to travel from the planet to our eyes. This is done only once, after the first approximation, but before the 2nd approximation. It will not be done between any other successive approximations.

For i=1 to 3

tc(i) = t(i) - A p(i)

Where A is the time (in days) required for light to travel 1 astronomical unit.

tau(1) = k {tc(3)-tc(2)}

tau(2) = k {tc(3)-tc(1)}

tau(3) = k {tc(2)-tc(1)}

nc(1) = tau(1)/tau(2)

nc(3) = tau(3)/tau(2)

**Recursive Procedure for Successive Approximations.**For i=1 to 3

x(i) = a(i) p(i) - Xs(i)

y(i) = b(i) p(i) - Ys(i)

z(i) = c(i) p(i) - Zs(i)

K(1) = { 2 [ r(2) r(3) + x(2) x(3) + y(2) y(3) + z(2) z(3) ] }^0.5

K(2) = { 2 [ r(1) r(3) + x(1) x(3) + y(1) y(3) + z(1) z(3) ] }^0.5

K(3) = { 2 [ r(1) r(2) + x(1) x(2) + y(1) y(2) + z(1) z(2) ] }^0.5

h(1) = tau(1)^2 / { K(1)^2 [ K(1)/3 + ( r(2)+r(3) )/2 ] }

h(2) = tau(2)^2 / { K(2)^2 [ K(2)/3 + ( r(1)+r(3) )/2 ] }

h(3) = tau(3)^2 / { K(3)^2 [ K(3)/3 + ( r(1)+r(2) )/2 ] }

For i=1 to 3

Begin

Queasy1 = (11/9) h(i)

Queasy2 = Queasy1

Repeat

Queasy3 = Queasy2

Queasy2 = Queasy1 / ( 1 + Queasy2)

Until abs( Queasy2 - Queasy3 ) / Queasy3 < 1E-12;

Yippy(i) = 1 + (10/11) Queasy2

End

n(1) = nc(1) yippy(2) / yippy(1)

n(3) = nc(3) yippy(2) / yippy(3)

nu(1) = nc(1) r(2)^3 { yippy(2) / yippy(1) - 1 }

nu(3) = nc(3) r(2)^3 { yippy(2) / yippy(3) - 1 }

K0 = { d(2) - d(1) nc(1) - d(3) nc(3) } / D

L0 = { d(1) nu(1) + d(3) nu(3) } / D

p(2) = K0 - L0 / r(2)^3

r(2) = { Re(2)^2 + q(2) p(2) + p(2)^2 }^0.5

f{r(2)} = r(2)^8 - { Re(2)^2 + q(2) K0 + K0^2 } r(2)^6 + L0 {q(2) + 2 K0} r(2)^3 - L0^2

f{r(2)} = 0

df/dr(2) = 8 r(2)^7 - 6 { Re(2)^2 + q(2) K0 + K0^2 } r(2)^5 + 3 L0 {q(2) + 2 K0} r(2)^2

The customary initial guess at r(2)...

r2(j=0) = 1.0

Repeat over j

r2(j+1) = r2(j) - f{r2(j)} / df/dr2(j)

Until abs{ r2(j+1) - r2(j) } < 1.0E-12

Assign to r(2) the converged value from the loop:

r(2) = r2(j+1)

p(2) = K0 - L0 / r(2)^3

Q1 = n(1) Xs(1) - Xs(2) + n(3) Xs(3) + a(2) p(2)

Q2 = n(1) Ys(1) - Ys(2) + n(3) Ys(3) + b(2) p(2)

Q3 = n(1) Zs(1) - Zs(2) + n(3) Zs(3) + c(2) p(2)

Q4 = a(1) - a(3) c(1) / c(3)

Q5 = b(1) - b(3) a(1) / a(3)

Q6 = b(3) - b(1) a(3) / a(1)

Q7 = a(3) - a(1) c(3) / c(1)

p(1) = { [ Q1 - a(3) Q3 / c(3) ] / [ n(1) Q4 ] + [ Q2 - b(3) Q1 / a(3) ] / [n(1) Q5] } / 2

p(3) = { [ Q2 - b(1) Q1 / a(1) ] / [ n(3) Q6 ] + [ Q1 - a(1) Q3 / c(1) ] / [n(3) Q7] } / 2

r(1) = { Re(1)^2 + q(1) p(1) + p(1)^2 }^0.5

r(3) = { Re(3)^2 + q(3) p(3) + p(3)^2 }^0.5

Repeat this whole section (Recursive Procedure for Successive Approximations) until your values for r(i) converge.

MORE ORBIT DETERMINATION IN POSTS TO FOLLOW.

**Example (after Dubyago with improved accuracy).****Initial data.**In practice, what you determine observationally about a planet's position is right ascension and declination. So I'll present those as initial data and calculate a(i), b(i), c(i) from them. Remember that Xs,Ys,Zs is the position of the sun, while RA and DEC pertain to the planet whose orbit we're trying to determine. (By the way, the planet is asteroid 1590 Tsiolkovskaja, also known as 1933 NA.)t(1) = JD 2427255.460417 = 1 July 1933, 23h 3m 0s UT

Xs(1) = -0.169709

Ys(1) = +0.919710

Zs(1) = +0.398865

RA(1) = 19.46730000 hours

DEC(1) = -13.86869444 degrees

t(2) = JD 2427283.391181 = 29 July 1933, 21h 23m 18s UT

Xs(2) = -0.600429

Ys(2) = +0.751016

Zs(2) = +0.325697

RA(2) = 19.06218056 hours

DEC(2) = -14.11902778 degrees

t(3) = JD 2427312.342083 = 27 August 1933, 20h 12m 36s UT

Xs(3) = -0.908371

Ys(3) = +0.405220

Zs(3) = +0.175716

RA(3) = 18.98696667 hours

DEC(3) = -15.24394444 degrees

**Unit vectors from observer to planet.**a(1) = +0.363835156

b(1) = -0.900093900

c(1) = -0.239697622

a(2) = +0.266215600

b(2) = -0.932536300

c(2) = -0.243937090

a(3) = +0.246531192

b(3) = -0.932786462

c(3) = -0.262929245

**First Approximation.**tau(1) = 0.498016281

tau(2) = 0.978484047

tau(3) = 0.480467766

nc(1) = 0.508967195

nc(3) = 0.491032805

nu(1) = 0.060177805

nu(3) = 0.059462580

d(1) = 0.015443444

d(2) = 0.018648221

d(3) = 0.017700437

D = 0.001964676

Re(1) = 1.016740339

Re(2) = 1.015193850

Re(3) = 1.010058035

q(1) = 1.970356907

q(2) = 1.879285653

q(3) = 1.296252781

K0 = 1.067106795

L0 = 1.008749516

The Lagrange equation is

f(r2) = r(2)^8 - 4.174733956 r(2)^6 + 4.048615418 r(2)^3 - 1.017575585

The first derivative of the Lagrange equation is

f'(r2) = 8 r(2)^7 - 25.048403737 r(2)^5 + 12.145846255 r(2)^2

This actually takes a while to converge with Newton's method, so I'd consider going to Danby's method for faster convergence. Anyway, it converges to

r(2) = 1.898670742

p(2) = 0.919728216

n(1) = 0.517759189

n(3) = 0.499720304

Q1 = +0.303475173

Q2 = -0.930010982

Q3 = -0.255727953

Q4 = +0.139086706

Q5 = +0.476529097

Q6 = -0.322891519

Q7 = -0.152567065

p(1) = 0.884503909

p(3) = 1.110851828

r(1) = 1.886503769

r(3) = 1.922018156

**Correction for planetary abberation.**tc(1) = JD 2427255.455309

tc(2) = JD 2427283.385869

tc(3) = JD 2427312.335667

tau(1) = 0.497997293

tau(2) = 0.978461559

tau(3) = 0.480464267

nc(1) = 0.508959486

nc(3) = 0.491040514

**Second Approximation.**x(1) = +0.491522618

y(1) = -1.715846574

z(1) = -0.610878483

x(2) = +0.845274999

y(2) = -1.608695948

z(2) = -0.550052825

x(3) = +1.182230625

y(3) = -1.441407547

z(3) = -0.467791433

K(1) = 3.801232986

K(2) = 3.732555561

K(3) = 3.766593197

h(1) = 0.005401695

h(2) = 0.021826229

h(3) = 0.005168609

Yippy(1) = 1.005962773

Yippy(2) = 1.023636797

Yippy(3) = 1.005707071

nu(1) = 0.061204833

nu(3) = 0.059919537

n(1) = 0.517901529

n(3) = 0.499794774

K0 = 1.067097940

L0 = 1.020939404

The Lagrange equation is

f(r2) = r(2)^8 - 4.174698414 r(2)^6 + 4.097521445 r(2)^3 - 1.042317267

The first derivative of the Lagrange equation is

f'(r2) = 8 r(2)^7 - 25.048190484 r(2)^5 + 12.292564335 r(2)^2

Newton's method converges to...

r(2) = 1.896390493

p(2) = 0.917399712

[Skipping the Q's.]

p(1) = 0.882742391

p(3) = 1.107232428

r(1) = 1.884757972

r(3) = 1.918706335

**Higher Approximations.**You've seen how it works. You just keep taking p(i) from the end of the latest approximation and plugging it back into the top of the next approximation and repeating the math with the improved numbers. You keep repeating the recursive section until all the values for r(i) have converged.

Jerry Abbott

Last edited: