Rational simplification bug ("Quotient by a polynomial of higher degree" error message)

I am using maxima-5.26 indirectly via sage-5.4.1, and I have only had
a week's worth of experience with sage (and indirectly maxima).
But it's been a productive week since I have been able to put together
a sage script that confirms certain symbolic sums over nbody indices (where
nbody is the number of gravitational bodies considered in a
General relativity calculation.)

The confirmation result should be zero and is stored in a Residual
Vector with just 3 (spatial) components.  Here is a typical good
result for nbody=2.

sage: #attach script (that is attached to this post) that does most of the work
sage: attach "location_check.sage"

sage: Residual[0]
(2*(beta + gamma)*G*M_1/(sqrt((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 -
x_1)^2)*c^2) - (vx_0^2 + vy_0^2 + vz_0^2)*gamma/c^2 - 1)*(x_0 -
x_1)*G*M_1/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2)^(3/2) +
1/2*(4*gamma + 3)*G*M_1*a_x1/(sqrt((z_0 - z_1)^2 + (y_0 - y_1)^2 +
(x_0 - x_1)^2)*c^2) + (x_0 - x_1)*G*M_1/((z_0 - z_1)^2 + (y_0 - y_1)^2
+ (x_0 - x_1)^2)^(3/2) - 1/2*(x_0 - x_1)*(2*(gamma + 1)*(vx_1^2 +
vy_1^2 + vz_1^2) - 4*(gamma + 1)*(vx_0*vx_1 + vy_0*vy_1 + vz_0*vz_1) -
2*(2*beta - 1)*U_1 - (z_0 - z_1)*a_z1 - (y_0 - y_1)*a_y1 - (x_0 -
x_1)*a_x1 - 3*((z_0 - z_1)*vz_1 + (y_0 - y_1)*vy_1 + (x_0 -
x_1)*vx_1)^2/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 -
x_1)^2))*G*M_1/(((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 -
x_1)^2)^(3/2)*c^2) + (vx_0 - vx_1)*(2*(gamma + 1)*((z_0 - z_1)*vz_0 +
(y_0 - y_1)*vy_0 + (x_0 - x_1)*vx_0) - (2*gamma + 1)*((z_0 - z_1)*vz_1
+ (y_0 - y_1)*vy_1 + (x_0 - x_1)*vx_1))*G*M_1/(((z_0 - z_1)^2 + (y_0 -
y_1)^2 + (x_0 - x_1)^2)^(3/2)*c^2) + 1/2*(4*(gamma + 1)*((z_0 -
z_1)*(vz_0 - vz_1) + (y_0 - y_1)*(vy_0 - vy_1) + (x_0 - x_1)*(vx_0 -
vx_1))*G*M_1*vx_1/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 -
x_1)^2)^(3/2) - 4*(gamma + 1)*G*M_1*a_x1/sqrt((z_0 - z_1)^2 + (y_0 -
y_1)^2 + (x_0 - x_1)^2) - (x_0 - x_1)*(2*(2*gamma + 1)*G*M_1/sqrt((z_0
- z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2) + vx_0^2 + vy_0^2 +
vz_0^2)*G*M_1/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2)^(3/2) -
2*((2*gamma + 1)*((z_0 - z_1)*(vz_0 - vz_1) + (y_0 - y_1)*(vy_0 -
vy_1) + (x_0 - x_1)*(vx_0 - vx_1))*G*M_1/((z_0 - z_1)^2 + (y_0 -
y_1)^2 + (x_0 - x_1)^2)^(3/2) + (z_0 - z_1)*G*M_1*vz_0/((z_0 - z_1)^2
+ (y_0 - y_1)^2 + (x_0 - x_1)^2)^(3/2) + (y_0 - y_1)*G*M_1*vy_0/((z_0
- z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2)^(3/2) + (x_0 -
x_1)*G*M_1*vx_0/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 -
x_1)^2)^(3/2))*vx_0)/c^2 - 1/2*(2*(2*beta - 1)*(x_0 -
x_1)*G^2*M_1^2/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2)^2 -
(2*gamma + 1)*(x_0 - x_1)*(vx_0^2 + vy_0^2 + vz_0^2)*G*M_1/((z_0 -
z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2)^(3/2) - (x_0 - x_1)*(2*(gamma
+ 1)*(vx_1^2 + vy_1^2 + vz_1^2) - 2*(2*beta - 1)*U_1 - (z_0 -
z_1)*a_z1 - (y_0 - y_1)*a_y1 - (x_0 - x_1)*a_x1 - ((z_0 - z_1)*vz_1 +
(y_0 - y_1)*vy_1 + (x_0 - x_1)*vx_1)^2/((z_0 - z_1)^2 + (y_0 - y_1)^2
+ (x_0 - x_1)^2))*G*M_1/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 -
x_1)^2)^(3/2) + (2*(x_0 - x_1)*((z_0 - z_1)*vz_1 + (y_0 - y_1)*vy_1 +
(x_0 - x_1)*vx_1)^2/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2)^2
- 2*((z_0 - z_1)*vz_1 + (y_0 - y_1)*vy_1 + (x_0 -
x_1)*vx_1)*vx_1/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2) -
a_x1)*G*M_1/sqrt((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2) +
4*(gamma + 1)*((x_0 - x_1)*G*M_1*vx_0*vx_1/((z_0 - z_1)^2 + (y_0 -
y_1)^2 + (x_0 - x_1)^2)^(3/2) + (x_0 - x_1)*G*M_1*vy_0*vy_1/((z_0 -
z_1)^2 + (y_0 - y_1)^2 + (x_0 - x_1)^2)^(3/2) + (x_0 -
x_1)*G*M_1*vz_0*vz_1/((z_0 - z_1)^2 + (y_0 - y_1)^2 + (x_0 -
x_1)^2)^(3/2)))/c^2

sage: Residual[0].simplify_rational()
0

and similarly for Residual[1] and Residual[2].  Those simplifications
take just a second or so on an entry level PC.

So all is well for
nbody=2.  When I try the same script except for changing nbody=2 to
nbody=3, then the result for Residual[0] is roughly twice as long
(involving additional variables like M2, x_2, y_2, z_2, etc.) and
therefore presumably more difficult to simplify to zero.  In fact,
I currently get the following maxima error message from that attempt:

sage: Residual[0].simplify_rational()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/irwin/<ipython console> in <module>()

/home/software/sage/sage-5.4.1/local/lib/python2.7/site-packages/sage/misc/decorators.pyc in wrapper(*args, **kwds)
     689                     kwds[new_name] = kwds[old_name]
     690                     del kwds[old_name]
--> 691             return func(*args, **kwds)
     692
     693         return wrapper

/home/software/sage/sage-5.4.1/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.simplify_rational (sage/symbolic/expression.cpp:30806)()

/home/software/sage/sage-5.4.1/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, x, name)
     193
     194         if isinstance(x, basestring):
--> 195             return cls(self, x, name=name)
     196         try:
     197             return self._coerce_from_special_method(x)

/home/software/sage/sage-5.4.1/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __init__(self, parent, value, is_name, name)
     620                 self._name = parent._create(value, name=name)
     621             except (TypeError, KeyboardInterrupt, RuntimeError, ValueError), x:
--> 622                 raise TypeError, x
     623
     624     def _latex_(self):

TypeError: ECL says: Error executing code in Maxima: Quotient by a polynomial of higher degree

I did a google search for that last error message, and there appeared
to be a number of hits for older maxima versions that were due to bugs
for those versions.  If anybody here thinks they know what the Maxima
problem is for this case of rational simplification, I have the
capability to rebuild sage with a patched version of maxima-5.26 to
try out your idea for a fix.  However, you might prefer to figure this
out directly.  Therefore, I have attached the relatively short sage
script to give you the opportunity to (1) verify all is well for the
rational simplification of Residual[0] generated by that script with
nbody=2. Such a check should only take you a second or so.  However,
if you edit that script to change one line from nbody=2 ==> nbody=3,
you should be able to confirm the above rational simplification error
(which should serve as a good test case for the fix).

Alan
__________________________
Alan W. Irwin

Astronomical research affiliation with Department of Physics and Astronomy,
University of Victoria (astrowww.phys.uvic.ca).

Programming affiliations with the FreeEOS equation-of-state
implementation for stellar interiors (freeeos.sf.net); the Time
Ephemerides project (timeephem.sf.net); PLplot scientific plotting
software package (plplot.sf.net); the libLASi project
(unifont.org/lasi); the Loads of Linux Links project (loll.sf.net);
and the Linux Brochure Project (lbproject.sf.net).
__________________________

Linux-powered Science
__________________________
-------------- next part --------------
# Verify location-dependent term in TCB-TCG (or TDB - TT) relation.
# See research note for notation: 
#
# TCB - TCG = (r_o - r_E)/(86400*c^2) dot Vector + 
#             (1/c^2) integral_T_0^TCB Integrand dt
#
# where (r_o - r_E) is the vector position of the observer relative to
# geocentre, Vector is the time ephemeris vector (with velocity units)
# the integral is the time ephemeris integral (with
# units of days) and the integrand is the time ephemeris
# integrand (with squared velocity units).  From the research note we have
#
# Vector = (v_E + P_E/c^2)
# Integrand = (U_ext + v_E^2/2 + Q_E/c^2)
#
# The 4-dimensional transformation between BCRS and GCRS for which the
# above expression for TCB-TCG is the time part is done at the
# velocity, v_E, of the GCRS relative to the BCRS.  Furthermore,
# Vector and Integrand are functions of r_E(t) = x_E(t), y_E(t), z_E(t), and
# v_E(t) where t is coordinate time, and these time varying functions can be
# determined from planetary ephemerides.  The 
# (r_o - r_E)/(86400*c^2) dot Vector term in the above expression for TCB-TCG
# is simply the first term in a Taylor series approximation for the integral for
# positions that differ (slightly) from the geocentre.

# This sage routine verifies that 
# the expression for Vector stated in Fukushima 1995 and also
# in my research note should be consistent with the following relationship
#
# derivative Vector(t, r_E(t), V_E(t)) with respect to t = 
# gradient Integrand(r_E, v_E) with respect to r_E
#
# which must be followed for the first term in
# a Taylor series expansion of Integral.  Therefore proof that this relationship
# holds proves that Vector as stated in Fukushima 1995 and my
# research note actually does have the correct form.

# The remaining sage definitions for P_E, U_ext, Q_E are taken
# directly from the research note, and the equation for the point-mass
# acceleration of the Earth (that is derivative V_E(t) with respect to t), is
# given by Equation 8.1 from
# ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/ExplSupplChap8.pdf.

G, c, beta, gamma = var("G,c,beta, gamma")

U_ext_summand = []

# na_E = Newtonian acceleration of the Earth.
na_E_summand = []
na_E_summand.append([])
na_E_summand.append([])
na_E_summand.append([])

# This is the curly V_ext vector from the research note.
curly_V_ext_summand = []
curly_V_ext_summand.append([])
curly_V_ext_summand.append([])
curly_V_ext_summand.append([])

# Delta_ext is the scalar quantity from the research note.
Delta_ext_summand = []

# Need explicit number of bodies.  Use small number rather than 11 (or
# 9000 if you were including asteroids) to make the symbolic
# processing easier.  But if we obtain a proof for nbody=4, say, then
# the proof is correct for arbitrary numbers of bodies.
nbody=2

# independent time variable, t.
t=var("t")
to_time = {}
from_time = {}

# Earth can be any body so long as it is less than nbody.
E = 0

# Note that r_E is geocentre location, while r_obs is the observer
# location.  We take the gradient of Integrand with respect to
# observer location, but then take the limit of that result for r_obs
# -> r_E.  Some quantities in the integrand depend on r_obs, and some
# depend on r_E, and we have to carefully distinguish the two in order
# to get an accurate linear term in the Taylor series for the
# integrand as a function of the change in observer location from
# geocentre.
r_E = var("x_%d, y_%d, z_%d" % (E, E, E))
r_obs = var("x_obs, y_obs, z_obs")
to_geocentre = {}
# dictionary to help convert from observer location to geocentre.
to_geocentre[r_obs[0]] = r_E[0]
to_geocentre[r_obs[1]] = r_E[1]
to_geocentre[r_obs[2]] = r_E[2]

# update argument dictionaries with Earth position components.  These
# conversions are only done for the Vector which is the limit as r_obs
# -> r_E for the partial derivative of the Integral with respect
# to r_obs so r_obs does not appear in these argument transformations.
r_E_t = (function("x_%d_t" % E, t), function("y_%d_t" % E, t), function("z_%d_t" % E, t))
to_time[r_E[0]] = r_E_t[0]
to_time[r_E[1]] = r_E_t[1]
to_time[r_E[2]] = r_E_t[2]
from_time[r_E_t[0]] = r_E[0]
from_time[r_E_t[1]] = r_E[1]
from_time[r_E_t[2]] = r_E[2]
v_E = var("vx_%d, vy_%d, vz_%d" % (E, E, E))

# update argument dictionaries with Earth velocity components
v_E_t = (r_E_t[0].diff(t), r_E_t[1].diff(t), r_E_t[2].diff(t))
to_time[v_E[0]] = v_E_t[0]
to_time[v_E[1]] = v_E_t[1]
to_time[v_E[2]] = v_E_t[2]
from_time[v_E_t[0]] = v_E[0]
from_time[v_E_t[1]] = v_E[1]
from_time[v_E_t[2]] = v_E[2]

# a_E = Earth acceleration from explanatory supplement 8.1 (AS-8.1)
a_E_summand = []
a_E_summand.append([])
a_E_summand.append([])
a_E_summand.append([])

for i in range(nbody):
    if i != E:
        # position vector of ith planet.
        r_i = var("x_%d, y_%d, z_%d" % (i, i, i))
        
        # update argument dictionaries with non-Earth position components
        r_i_t = (function("x_%d_t" % i, t), function("y_%d_t" % i, t), function("z_%d_t" % i, t))
        to_time[r_i[0]] = r_i_t[0]
        to_time[r_i[1]] = r_i_t[1]
        to_time[r_i[2]] = r_i_t[2]
        from_time[r_i_t[0]] = r_i[0]
        from_time[r_i_t[1]] = r_i[1]
        from_time[r_i_t[2]] = r_i[2]

        # IMPORTANT. Use observer location rather than geocentre for
        # all terms here.  This effectively means all terms in the
        # time ephemeris Integrand other than the planetary
        # accelerations and potentials which are treated symbolically
        # here in any case.  For the time ephemeris Vector, the
        # transformation from observer location -> geocentre is taken
        # in any case for all quantities.

        # position vector of the observer relative to the ith planet.
        r_obs_i = [r_obs[0]-r_i[0], r_obs[1]-r_i[1], r_obs[2]-r_i[2]]
        # distance between observer and ith planet.
        dist_obs_i = sqrt(r_obs_i[0]^2 + r_obs_i[1]^2 + r_obs_i[2]^2)

        mu_i = G*var("M_%d" % i)

        U_ext_summand.append(mu_i/dist_obs_i)
        
        na_E_summand[0].append(-mu_i*r_obs_i[0]/dist_obs_i^3)
        na_E_summand[1].append(-mu_i*r_obs_i[1]/dist_obs_i^3)
        na_E_summand[2].append(-mu_i*r_obs_i[2]/dist_obs_i^3)
    
        # velocity of the ith planet
        v_i = var("vx_%d, vy_%d, vz_%d" % (i, i, i))

        # update argument dictionaries with non-Earth velocity components
        v_i_t = (r_i_t[0].diff(t), r_i_t[1].diff(t), r_i_t[2].diff(t))
        to_time[v_i[0]] = v_i_t[0]
        to_time[v_i[1]] = v_i_t[1]
        to_time[v_i[2]] = v_i_t[2]
        from_time[v_i_t[0]] = v_i[0]
        from_time[v_i_t[1]] = v_i[1]
        from_time[v_i_t[2]] = v_i[2]

        curly_V_ext_summand[0].append(mu_i*v_i[0]/dist_obs_i)
        curly_V_ext_summand[1].append(mu_i*v_i[1]/dist_obs_i)
        curly_V_ext_summand[2].append(mu_i*v_i[2]/dist_obs_i)
        
        # U_i and a_i are the symbolic gravitational potential and 
        # acceleration observed at the position of the ith planet.
        U_i = var("U_%d" % i)
        a_i = []
        a_i.append(var("a_x%d" % i))
        a_i.append(var("a_y%d" % i))
        a_i.append(var("a_z%d" % i))
        
        # update from_time dictionary with non-Earth acceleration components
        a_i_t = (r_i_t[0].diff(t, 2), r_i_t[1].diff(t, 2), r_i_t[2].diff(t, 2))
        from_time[a_i_t[0]] = a_i[0]
        from_time[a_i_t[1]] = a_i[1]
        from_time[a_i_t[2]] = a_i[2]

        # Calculate each of the factors entering into the relationship
        # for Delta_ext in turn.
        r_obs_i_dot_v_i = r_obs_i[0]*v_i[0] + r_obs_i[1]*v_i[1] + r_obs_i[2]*v_i[2]
        r_obs_i_dot_a_i = r_obs_i[0]*a_i[0] + r_obs_i[1]*a_i[1] + r_obs_i[2]*a_i[2] 
        v_i_dot_v_i = v_i[0]*v_i[0] + v_i[1]*v_i[1] + v_i[2]*v_i[2]
        Delta_0 = r_obs_i_dot_v_i^2/(2*dist_obs_i^2)
        Delta_1 = r_obs_i_dot_a_i/2
        Delta_2 = -(1+gamma)*v_i_dot_v_i
        Delta_3 = -(1 - 2*beta)*U_i
        Delta_ext_summand.append(mu_i*(Delta_0 + Delta_1 + Delta_2 + Delta_3)/dist_obs_i)
        # Calculate factors entering into equation 7 of the research
        # note excluding the first two and the fourth (which are
        # calculated outside this loop).
        v_E_dot_v_i = v_E[0]*v_i[0] + v_E[1]*v_i[1] + v_E[2]*v_i[2]
        a_E_2 = (1 - 2*beta)*U_i
        a_E_4 = (1 + gamma)*v_i_dot_v_i
        a_E_5 = -2*(1 + gamma)*v_E_dot_v_i
        a_E_6 = -(3/2)*(r_obs_i_dot_v_i/dist_obs_i)^2
        a_E_7 = -(1/2)*r_obs_i_dot_a_i
        a_E_bracket_0 = (mu_i/(dist_obs_i^3*c^2))*(a_E_2+a_E_4+a_E_5+a_E_6+a_E_7)
        
        v_E_dot_r_obs_i = v_E[0]*r_obs_i[0] + v_E[1]*r_obs_i[1] + v_E[2]*r_obs_i[2]
        a_E_8 = 2*(1 + gamma)*v_E_dot_r_obs_i
        a_E_9 = -(1 + 2*gamma)*r_obs_i_dot_v_i
        a_E_bracket_1 = (mu_i/(dist_obs_i^3*c^2))*(a_E_8+a_E_9)
        a_E_bracket_2 = ((3 + 4*gamma)/2)*(mu_i/(dist_obs_i*c^2))
        
        a_E_summand[0].append(-a_E_bracket_0*r_obs_i[0] + \
                                   a_E_bracket_1*(v_E[0]-v_i[0]) + \
                                   a_E_bracket_2*a_i[0])
        a_E_summand[1].append(-a_E_bracket_0*r_obs_i[1] + \
                                   a_E_bracket_1*(v_E[1]-v_i[1]) + \
                                   a_E_bracket_2*a_i[1])
        a_E_summand[2].append(-a_E_bracket_0*r_obs_i[2] + \
                                   a_E_bracket_1*(v_E[2]-v_i[2]) + \
                                   a_E_bracket_2*a_i[2])
                                          
U_ext = sum(U_ext_summand)

na_E = []
na_E.append(sum(na_E_summand[0]))
na_E.append(sum(na_E_summand[1]))
na_E.append(sum(na_E_summand[2]))

# update from_time dictionary with Earth acceleration components
a_E_t = (r_E_t[0].diff(t, 2), r_E_t[1].diff(t, 2), r_E_t[2].diff(t, 2))
from_time[a_E_t[0]] = na_E[0].subs(to_geocentre)
from_time[a_E_t[1]] = na_E[1].subs(to_geocentre)
from_time[a_E_t[2]] = na_E[2].subs(to_geocentre)

curly_V_ext = []
curly_V_ext.append(sum(curly_V_ext_summand[0]))
curly_V_ext.append(sum(curly_V_ext_summand[1]))
curly_V_ext.append(sum(curly_V_ext_summand[2]))

Delta_ext = sum(Delta_ext_summand)

v_E_dot_v_E = v_E[0]*v_E[0] + v_E[1]*v_E[1] + v_E[2]*v_E[2]
v_E_dot_curly_V_ext = v_E[0]*curly_V_ext[0] + v_E[1]*curly_V_ext[1] + v_E[2]*curly_V_ext[2]

Q_E = ((1 - 2*beta)/2)*U_ext*U_ext + ((1 + 2*gamma)/2)*U_ext*v_E_dot_v_E + \
    v_E_dot_v_E^2/8 -2*(1 + gamma)*v_E_dot_curly_V_ext - Delta_ext
Integrand = U_ext + v_E_dot_v_E/2 + Q_E/(c*c)

Gradient = (Integrand.diff(r_obs[0]).subs(to_geocentre), Integrand.diff(r_obs[1]).subs(to_geocentre), Integrand.diff(r_obs[2]).subs(to_geocentre))

# Vector quantities
# Note, must keep Newtonian and Post-Newtonian parts of Vector
# separate because acceleration treated differently in the two
# cases.

# Post_Newtonian:
#P_E_mult_0 = (2 + gamma)*U_ext + v_E_dot_v_E/2
P_E_mult_0 = (1+2*gamma)*U_ext + v_E_dot_v_E/2
P_E_mult_1 = -2*(1 + gamma)
P_E = []
P_E.append(P_E_mult_0*v_E[0] + P_E_mult_1*curly_V_ext[0])
P_E.append(P_E_mult_0*v_E[1] + P_E_mult_1*curly_V_ext[1])
P_E.append(P_E_mult_0*v_E[2] + P_E_mult_1*curly_V_ext[2])

# Convert to geocentre, convert to time representation, differentiate
# wrt time, then convert that result back to normal representation.
P_E[0] = P_E[0].subs(to_geocentre).subs(to_time).diff(t).subs(from_time)
P_E[1] = P_E[1].subs(to_geocentre).subs(to_time).diff(t).subs(from_time)
P_E[2] = P_E[2].subs(to_geocentre).subs(to_time).diff(t).subs(from_time)

# Newtonian part of result is just v_E and its time derivative is
# a_E.

# Accumulate components of a_E including hitherto missing first,
# second, and fourth terms.
factor = 1 - (2*(beta+gamma)/c^2)*U_ext + (gamma/c^2)*v_E_dot_v_E
a_E = []
a_E.append(na_E[0]*factor + sum(a_E_summand[0]))
a_E.append(na_E[1]*factor + sum(a_E_summand[1]))
a_E.append(na_E[2]*factor + sum(a_E_summand[2]))

a_E[0] = a_E[0].subs(to_geocentre)
a_E[1] = a_E[1].subs(to_geocentre)
a_E[2] = a_E[2].subs(to_geocentre)

Vector = []
Vector.append(a_E[0] + P_E[0]/c^2)
Vector.append(a_E[1] + P_E[1]/c^2)
Vector.append(a_E[2] + P_E[2]/c^2)

Residual = (Vector[0]-Gradient[0], Vector[1]-Gradient[1], Vector[1]-Gradient[2])