I will go through 2 slightly different calculations of the problem below.
Jupiter emits roughly twice as much energy into space as it receives from the Sun. The difference is thought to be produced by continuing slow gravitational contraction of the planet. Assuming that the gravitational potential energy of the planet ' $GM^2/R$ (where the missing constant of proportionality is of order unity) estimate by how much the radius of the planet must change each year.
$$R_{\text{Jup}}=7 \times 10^7\,\mathrm{m}$$
$$d_{\text{Sun-Jup}}=5.2 \,\mathrm{AU}=5.2 \times 1.5 \times 10^{11}\,\mathrm{m}$$
$$L_{\text{Sun}}=4 \times 10^{26}\,\mathrm{W}$$
$$M_{\text{Jup}}=1.9 \times 10^{27}\,\mathrm{kg}$$
The Suns power, $L_{\text{Sun}}$ is spread over a sphere of radius $d_{\text{Sun-Jup}}$ given by $$F_{\text{Sun}}=\frac{L_{\text{Sun}}}{4\pi \,{d_{\text{Sun-Jup}}}^2}$$ where $F_{\text{Sun}}$ is the flux emitted by the sun. But Jupiter only intercepts $\pi R_{\text{Jup}}^2$ (a disk) of this flux, so the intercepted power (or luminosity received from the Sun) is
$$\pi R_{\text{Jup}}^2\frac{L_{Sun}}{4\pi \,{d_{\text{Sun-Jup}}}^2}=\frac{1}{4}L_{\text{Sun}}\left(\frac{R_{\text{Jup}}}{d_{\text{Sun-Jup}}}\right)^2=\frac{4\times 10^{26}}{4}\times\left(\frac{7\times 10^7}{5.2\times 1.5\times 10^{11}}\right)^2\approx 8\times 10^{17}\,\mathrm{W}$$
The excess power $P$ released through gravitational contraction is
$$P=\frac{d}{dt}\frac{GM^2}{R}=GM^2\frac{d}{dR}\frac{1}{R}\frac{dR}{dt}=-\frac{GM^2}{R^2}\frac{dR}{dt}$$ I wish to solve the differential equation $$P=-\frac{GM^2}{R^2}\frac{dR}{dt}\tag{1}$$
Applying boundary conditions such that $$t=\begin{cases} 0, &\text{when}\,\,\, R=7\times 10^7\,\mathrm{m}\,\,\, \text{(Jupiters initial radius)} \\ 3.1\times 10^7, &\text{when}\,\,\, R= \text{Jupiter radius after 1 year}\,\,(3.1\times 10^7\,\mathrm{s})\end{cases}$$
Separating variables and integrating $(1)$
$$\int_{R'=R_{Jup}}^{R}\frac{1}{{R'}^2}dR'=-\frac{P}{GM^2}\int_{t=0}^{1\,\mathrm{Yr}}dt$$
$$\implies-\bigg[\frac{1}{R'}\bigg]_{7\times 10^7}^{R}=-\frac{P}{GM^2}\bigg[t\bigg]_{0}^{3.1 \times 10^7}$$
$$\implies\frac{1}{R}-\frac{1}{7\times 10^7}=\frac{8\times 10^{17}\times 3.1\times 10^7}{6.7\times 10^{-11}\times\left(1.9\times 10^{27}\right)^2}$$
$$\implies\frac{1}{R}=\frac{8\times 10^{17}\times 3.1\times 10^7}{6.7\times 10^{-11}\times\left(1.9\times 10^{27}\right)^2}+\frac{1}{7\times 10^7}=\frac{7\times 10^7\times 8\times 10^{17}\times 3.1\times 10^7+6.7\times 10^{-11}\times\left(1.9\times 10^{27}\right)^2}{7\times 10^7\times 6.7\times 10^{-11}\times\left(1.9\times 10^{27}\right)^2}$$
$$\implies \frac{1}{R}=\frac{56\times 3.1\times 10^{31}+6.7\times 3.61\times 10^{43}}{7\times 6.7\times 3.61\times 10^{50}}=\frac{173.6+24.187\times 10^{12}}{169.309\times 10^{19}}$$
$$\implies R=\frac{169.309\times 10^{19}}{173.6+24.187\times 10^{12}}$$
Using Wolfram Alpha to evaluate this
this corresponds to a radius decrease (per year) of
$0.0005024186546455856502068684140705373813750095\,\mathrm{m}$
There is another way to solve this problem, and that is by considering the fractional change in $R$
Let $\dot R =\frac{dR}{dt}$ and rearrange $(1)$ such that
$$\frac{\dot{R}}{R}=-P\frac{R}{GM^2}$$
Inserting the same values as before to exactly the same precision
$$\frac{\dot{R}}{R}=-8\times 10^{17}\times\frac{7\times 10^7}{6.7\times 10^{-11}\times \left(1.9\times 10^{27} \right)^2}=-\frac{56\times 10^{24}}{6.7\times 3.61\times 10^{43}}$$ $$\implies \frac{\dot R}{R}=-\frac{56}{24.187\times 10^{19}}\tag{2}$$ $$\implies \frac{\dot R}{R}=-\frac{56}{24.187\times 10^{19}}\approx -2.3 \times 10^{-19}\,\mathrm{s}^{-1}$$
With $3.1 \times 10^{7}\,\mathrm{s}$ per year, the expected fractional contraction per year is about $7 \times 10^{−12}$.
Now, finally, to get the actual radius this corresponds to I multiply the RHS of $(2)$ by the number of seconds in a year and the initial radius of Jupiter
$$\implies R=\frac{56}{24.187\times 10^{19}}\times 7\times 10^7 \times 3.1\times 10^7$$
On simplification this gives
$$R=\frac{1215.2\times 10^{14}}{24.187\times 10^{19}}=\frac{1215.2}{24.187\times 10^5}$$
Evaluating this with Wolfram Alpha again:
So the radius decrease this time is $0.0005024186546491917145574068714598751395377682\,\mathrm{m}$
Close inspection shows that these two numbers (to $46$ decimal places) are not identical:
$R=0.000\color{blue}{50241865464}55856502068684140705373813750095\,\mathrm{m}$ for the first method
$R=0.000\color{blue}{50241865464}91917145574068714598751395377682\,\mathrm{m}$ for the second method
I used the same precision throughout and only computed with Wolfram Alpha at the end, so why are these numbers not exactly equal?
Or, put in another way,
Why is
$$\frac{169.309\times 10^{19}}{173.6+24.187\times 10^{12}}\ne\frac{1215.2}{24.187\times 10^5}?$$
Remark: I'm sorry for asking such a strange question. I just can't figure out why they are not equal and it's driving me mad. I also apologize for including many steps in the calculations, but I wanted to show each step to hopefully account for any possible errors.
In the first case you solve a differential equation $$ \dot R = -cR^2\implies \frac1{R(t_0+\Delta t)}-\frac1{R(t_0)}=cΔt \\ \frac{R(t_0+\Delta t)-R(t_0)}{R(t_0)}=-\frac{cR(t_0)Δt}{1+cR(t_0)Δt} =-cR(t_0)Δt+(cR(t_0)Δt)^2-(cR(t_0)Δt)^3\pm\dots $$ In the second case you use the derivative value at $t_0$ to compute the approximation $R_1\approx R(t_1)$ at $t_1=t_0+Δt$, which in effect is an Euler step $$ R_1= R_0+\dot R(t_0)Δt\implies\frac{R_1-R_0}{R_0}= -cR_0Δt $$ This is in first order the same, but misses the higher order terms. As you computed, $cR_0Δt\approx 7\cdot 10^{-12}$, so that the first 11 non-zero digits of both results should be the same, differences are to be expected at least in the 13th digit.