Q:

Accurate Equation for Earth's Gravitational Binding Energy

This is a relatively important question for anyone who can answer it. I am trying to find the equation that accurately solves for Earth's Gravitational Binding Energy. The information below is from the wikipedia page:
Assuming that the Earth is a uniform sphere (which is not correct, but is close enough to get an order-of-magnitude estimate) with M = 5.97 x 10^24 kg and r = 6.37 x 10^6 m, U is 2.24 x 10^32 J. This is roughly equal to one week of the Sun's total energy output. It is 37.5 MJ/kg, 60% of the absolute value of the potential energy per kilogram at the surface.
The actual depth-dependence of density, inferred from seismic travel times (see Adams–Williamson equation), is given in the Preliminary Reference Earth Model (PREM). Using this, the real gravitational binding energy of Earth can be calculated numerically as U = 2.487 x 10^32 J.
So what I wish to not for the latter results (2.487e+32 J) is what is the actual equation is used to get this result. 
And I do not mean the standard GBE equation which gave the other result above.

A:

I am trying to find the equation that accurately solves for Earth's Gravitational Binding Energy.

There isn't a single equation for this (and much of real science does not yield convenient single formulas as solutions).  What the PREM produces for density is a set of piecewise approximate polynomial functions that model the theorized density that produces a reasonable match to measured data (like for example seismic data).
So working from this page by Dave Typinski you could develop an integral equation for the gravitational binding energy.  I'm not going to actually do that myself, but the density polynomial functions are just quadratic, and you can apply this simply enough to the integral for gravitational binding energy :
$$U = 16\pi^2 G \int_0^R r \rho(r) \left[ \int_0^r \rho(r) r^2 dr \right] dr $$
This is tedious to do with the piecewise functions for density, but not difficult.
For a constant density $\rho(r)=\rho_0$ you can see this reduces to the familiar equation :
$$U = 16\pi^2G \rho_0^2 \int_0^R r \left[\frac 1 3 r^3\right] dr = \frac {16\pi^2G} {15} \rho_0^2 R^5 = \frac {3}{5} \frac {GM^2} R$$

So what I wish to not for the latter results ($2.487\times 10^{32}\,J$) is what is the actual equation is used to get this result. 

I can't answer that definitively, but it would would most likely be either a numerical integration or something like I've outlined.

A:

To calculate the gravitational energy should use two integrals
$M(r)=4\pi\int_0^r\rho(r')r'^2dr'$
$U=4\pi G\int_0^{R_E}\rho (r)M(r)rdr$
$R_E=6.371×10^6 m$ 
We use the density distribution obtained on the basis of PREM (see Earth's Gravitational Binding Energy) 

As a result, we find $U=2.48848×10^{32} J$ which is consistent with the value $2.487×10^{32} J$ shown on the wikipedia page  and calculated using PREM. 
Below is the Wolfram code debugged for Mathematica 11.3.
A = {{0, 0, 0, 0}, {1.2215*10^6, -2.1773*10^-10, 1.9110*10^-8, 
1.3088*10^4}, {3.4800*10^6, -2.4123*10^-10, 1.3976*10^-4, 
1.2346*10^4}, {3.6300*10^6, 0, -5.0007*10^-4, 
7.3067*10^3}, {5.7010*10^6, -3.0922*10^-11, -2.4441*10^-4, 
6.7823*10^3}, {5.7710*10^6, 0, -2.3286*10^-4, 
5.3197*10^3}, {5.9710*10^6, 0, -1.2603*10^-3, 
1.1249*10^4}, {6.1510*10^6, 0, -5.9706*10^-4, 
7.1083*10^3}, {6.3466*10^6, 0, 1.0869*10^-4, 
2.6910*10^3}, {6.3560*10^6, 0, 0, 2.9000*10^3}, {6.3680*10^6, 0, 
0, 2.6000*10^3}, {6.3710*10^6, 0, 0, 1.0200*10^3}};

rho[x_] := 
 Piecewise[
  Table[{A[[i, 2]]*x^2 + A[[i, 3]]*x + A[[i, 4]], 
    A[[i - 1, 1]] < x <= A[[i, 1]]}, {i, 2, Length[A]}]]    
    Plot[rho[x], {x, 0, 6.3710*10^6}, PlotRange -> All, PlotPoints -> 100,
      AxesOrigin -> {0, 0}, 
     AxesLabel -> {"r, m", 
       "\[Rho], kg/\!\(\*SuperscriptBox[\(m\), \(3\)]\)"}]

G = 6.67408*10^\[Minus]11; R = 6.3710*10^6;
4*Pi*NIntegrate[
  rho[x]*x^2, {x, 0, 
   R}](*a total Earth mass; the NASA figure of 5.9736\[Times]10^24 \
kg.*)

M = 
  Interpolation[
   Table[{r, 4*Pi*NIntegrate[rho[x]*x^2, {x, 0, r}]}, {r, 1, R, 
     10^4}]];
Plot[M[x], {x, 0, R}, PlotRange -> All, AxesOrigin -> {0, 0}, 
 AxesLabel -> {"r, m", "M, kg"}]

 U = 
 4*Pi*G*NIntegrate[
   rho[r]*r*M[r], {r, 0, R}](*Earth's Gravitational Binding Energy*)

(*Out[]= 2.48848*10^32*)

