I'm not sure I understand your code or what (1) and (2) above means.

So, what is the automatic default rule for this integration? From the manual, there's "Trapezoidal Rule", "Multipanel Rule" and so on. How do I know which rule that Mathematica chose?

Sometimes there's a way to tell which numerical method Mathematica used by using a trace of some kind. Usually there isn't. The method chosen might be a mixture of multiple methods or could even be a method not documented. It could be anything. If you don't specify a value for the "Method" option, the function is free to do whatever it wants to try and solve your problem. It doesn't try to communicate what it did to solve it - often what it did isn't straightforward enough anyway.

If you want to try to examine the output of the function, you can sometimes determine what numerical method was used, although I don't recommend it. http://mathematica.stackexchange.com/questions/145/how-to-find-out-which-method-mathematica-selected

Integrate integrates the polynomials symbolically and evaluates them at the endpoints of their subintervals; the result of integrating a standard order-3 InterpolatingFunction will be an order-4 InterpolatingFunction, with the degree-4 polynomials being the integrals of the degree-3 polynomials satisfying the appropriate initial conditions. It does not use a numerical method at all.

No problem.

Here are the steps I have in mind:

1. NDSolve produces y -> InterpolatingFunction[..]. The InterpolatingFunction will be order 3. Order 3 means the interpolation consists of piecewise degree-3 polynomials that interpolate between the steps of the solution computed by NDSolve.

2. If you integrate the solution with intsol = Integrate[y[x] /. First[sol], x], where sol = NDSolve[..] is the solution, then intsol, the value of the result of Integrate, will be of the form InterpolatingFunction[..][x], where the interpolating function is an order-4 interpolation, which consists of piecewise degree-4 polynomials that are the integrals of the degree-3 polynomials in the solution.

There is no second application of Integrate that I have in mind. If you want to use Integrate twice, then the order increases by one at each step.

By "the appropriate initial conditions", I mean on each subinterval ${(x_{j-1}, x_j)}$, Integrate performs the integral

$$I_{j-1} + \int_{x_{j-i}}^x p_j(t) \; dt$$

where the initial value $I_{j-1}$ is the accumulation the integrals over the preceding intervals..

If I understand what you did, Integrate, the symbolic integrator, was applied to the InterpolatingFunction. The interpolating function is basically a piecewise polynomial function (defined by data computed by NDSolve). Integrate effectively integrates each polynomial using the power rule and the rules of definite integrals, returning a function f such that f[x] equals the integral of the solution from the beginning of the domain up to x. It does not use numerical techniques.

Thank you Michael Rogers, this have given me an idea of how the integration works.

It's just that, for integrations without using Mathematica, we would use Simpson's Rules or Trapeziodal rule for integration for numerical integration, to obtain the value for the initial value before it. But in Mathematica without numerical integration, the InterpolatingFunction only show a graph or a plot (without any polynomials shown). Well I did obtain an accurate value for the integration of the InterpolatingFunction.

There's the set of rules in Mathematica for integration. So, now I need to know which is the correct integration rule that is default for Mathematica.

tutorial/NIntegrateIntegrationRules

I believe Mathematica automatic chooses the rule for the Integration, so what command do I need to type to find out which method it used?