Symbolic integration fails while numerical integration succeeds

Believe the numerical one. Mathematica simply could not do the symbolic integration. Symbolic integration will travel via a different code path. Here the symbolic integration done using Maple, and it agrees with the numerical solution given by Mathematica's NIntegrate

The analytical answer is

 (7/18)*hypergeom([-1/3, 1, 1], [2, 2], 7/8)-(7/36)*hypergeom([1/3, 1, 1], [2, 2], 7/8)

You can do it in Mathematica like this:

z = Integrate[((r^3 - 7)^(2/3)*(1 - (r^3 - 7)^(2/3)/r^2))/r^3, r];
r1 = Limit[z, r -> Infinity]
(* 3/4 + Pi/(3*Sqrt[3]) *)
r2 = Limit[z, r -> 2]
(r1 - r2) // N
Out[6]= 0.139837 + 0. I

Update:

There is some strange kernel buffering issue here. The integral will work, but it depends on order of thing. On a new FRESH kernel, I found the following: If I run the integral once (but the indefinite version), and then run it again, but use the definite integral now, it work ! Here is screen shot

I think what happens is that some definitions get loaded first time which did not exist before. Once these are loaded, then the next time it works.

You can do the analytic integral in Mathematica too, by telling it to perform the upper integration limit as follows:

With[
 {
  i = 
   Integrate[((r^3 - 7)^(2/3)*(1 - (r^3 - 7)^(2/3)/r^2))/r^3, r]
 },
 Simplify[
   Limit[i, r -> Infinity] - i /. r -> 2
  ]
 ]

(*
==> 23/64 + Pi/(3 Sqrt[3]) - 
 2 (-(1/7))^(1/3) Hypergeometric2F1[1/3, 1/3, 4/3, 8/7] + 
 2 (-(1/7))^(2/3) Hypergeometric2F1[2/3, 2/3, 5/3, 8/7]
*)

N[%]//Chop

(* ==> 0.139837 *)

I hope folks don't mind, but I'd like to revisit this question. Even in 11.3, Mathematica's approach to the symbolic integration is not optimal, as it gives an unnecessarily messy result. The integrand is real-valued for $r \ge \sqrt[3]{7}$ , so Mathematica's output for the antiderivative from 2 to $\infty$ should be exactly so as well:

integrand = ((r^3 - 7)^(2/3)*(1 - (r^3 - 7)^(2/3)/r^2))/r^3 //FullSimplify;
Plot[integrand, {r, 0, 20}, PlotRange -> {{0, 20}, {0, .12}}] 

Instead, as can be seen from the prior answers, Mathematica adds a complex integration constant that cancels out when the definite integral is evaluated numerically, leaving an imaginary term that is indistinguishable from zero, within the numerical precision of the operation; but that imaginary term shouldn't be there at all. A significantly cleaner result can be obtained with the Rubi integration package for Mathematica (http://www.apmaths.uwo.ca/~arich/):

rubi = Int[integrand, r]

And thus (Rubi doesn't do definite integrals, so the upper and lower bounds must be determined individually):

rubiU = Limit[rubi, r -> Infinity]
rubiL = rubi /. r -> 2
rubiU - rubiL
% // N

Essentially, compared to the Rubi result, Mathematica is adding an integration constant of $-0.3023-2.0944 i$:

mma = Integrate[integrand, r];
mma /. r -> 2.
rubi /. r -> 2.
%% - %

While the integration constant can have any constant value, this doesn't seem to be good program behavior; I don't see how it's any different from evaluating $\int x^2 \, dx$ and getting, say, $\frac{x^3}{3}-27+\frac{14 i}{7}$ rather than just $\frac{x^3}{3}$.