How to speed up the plot of NIntegrate?

You can express your integral in terms of a differential equation and use NDSolve. Since NDSolve builds up the solution as it goes, this is typically much faster.

Clear[y];
y[x_] = y[x] /. First[
  NDSolve[{y'[x] == Sin[x], y[0] == 0}, y[x], {x, 0, 10}]
];
t = AbsoluteTime[];
Plot[y[t], {t, 0, 10}]
AbsoluteTime[] - t

An alternative approach is to form an approximate interpolating function from your actual function, which is (often) cheaper to integrate. To wit,

nsin = FunctionInterpolation[Sin[x], {x, 0, 10}];
ni = Derivative[-1][nsin];

Plot[ni[t], {t, 0, 10}]

In this case, we know that the integral of sine can be expressed simply, so we can compare the exact function with the approximate one:

Plot[ni[t] - 2 Haversine[t], {t, 0, 10}, Frame -> True, PlotRange -> All]

Pretty darn good, if I do say so myself! If you need a bit more accuracy, FunctionInterpolation[] has a number of options you can tweak as needed; see the docs for details.

Can NIntegrate remember or make full use of the result of a smaller upper-limit integral? Or generally, how to speed up the plot involving NIntegrate or is there any principle to do it?

Below is a solution in the spirit of this request.

In short, we find adaptively sampled points, compute integral estimates over the intervals having those points as boundaries, compute cumulative sums, and list plot.

General procedure

First, note that NIntegrate is similar to Plot -- both use sampling points adaptively to produce their output. (E.g. see Max(Plot)Points and MaxRecursion.)

1. Use NIntegrate over the whole range the integral values to be plotted within. Gather the sampling points.

2. Partitition the sampling points into disjoint or overlapping intervals. Sort them accordingly to their boundaries.

3. For each interval compute an integral estimate using an integration rule. (Preserving the order of the intervals.)

4. Compute the cumulative sums of the integral estimates and pair them with the interval boundaries.

5. (List)Plot the obtained pairs.

Note that because of step 1 we will have a reasonably adaptive behavior.

How to obtain integration rules and integral estimates with them is shown in NIntegrate's advanced documentation. (It is also given below.)

Implementation

This implementation also allows sampling points specification instead of using NIntegrate's sampling points in step 1.

{absc, weights, errweights} = 
  NIntegrate`GaussKronrodRuleData[5, MachinePrecision];

IRuleEstimate[f_, {a_, b_}] := 
 Module[{integral, 
   error}, {integral, 
    error} = (b - a) Total@
     MapThread[{f[#1] #2, f[#1] #3} &, {Rescale[absc, {0, 1}, {a, b}], 
       weights, errweights}];
  {integral, Abs[error]}]

Clear[PPartialIntegrations]
Options[PPartialIntegrations] = {MinRecursion -> 1, PrecisionGoal -> 3, "PointsPerInterval" -> 2};
PPartialIntegrations[f_, {a_?NumericQ, b_?NumericQ}, 
   opts : OptionsPattern[]] :=
  Block[{minrec, pgoal, npint, res, x},
   minrec = OptionValue[MinRecursion];
   pgoal = OptionValue[PrecisionGoal];
   If[TrueQ[pgoal === Automatic], pgoal = 3];
   npint = OptionValue["PointsPerInterval"];
   If[npint === Automatic, npint = 2];
   res = Reap@
     NIntegrate[f[x], {x, a, b}, MinRecursion -> minrec, 
      PrecisionGoal -> pgoal, EvaluationMonitor :> Sow[x]];
   PPartialIntegrations[f, {a, b}, Sort[res[[2, 1]]][[1 ;; -1 ;; npint - 1]]]
   ];

PPartialIntegrations[f_, {a_?NumericQ, b_?NumericQ}, stepArg_] :=
  Block[{step = stepArg},
    If[step === Automatic, step = (b - a)/200];
    PPartialIntegrations[f, {a, b}, Range[a, b, step]]
    ] /; AtomQ[stepArg];

PPartialIntegrations[f_, {a_?NumericQ, b_?NumericQ}, 
   range : {_?NumericQ ..}] :=
  Block[{intervals, integrals, cumSums},
   res = Map[{#, IRuleEstimate[f, #][[1]]} &, Partition[range, 2, 1]];
   intervals = res[[All, 1]];
   integrals = res[[All, 2]];
   cumSums = Prepend[Accumulate[integrals], 0];
   Transpose[{Prepend[intervals[[All, 2]], intervals[[1, 1]]], cumSums}]
   ];

Basic examples

We can see that the implementation is reasonably fast.

AbsoluteTiming[PPartialIntegrations[Sin, {0, 10}, Automatic];]
(* {0.009395, Null} *)

AbsoluteTiming[
 PPartialIntegrations[Sin, {0, 10}, MinRecursion -> 2, PrecisionGoal -> 5];]
(* {0.006568, Null} *)

ListLinePlot[PPartialIntegrations[Sin, {0, 10}, Automatic]]

Advanced example

Consider the function:

Plot[Sin[x] + 2 Exp[-100 (x - 5)^2], {x, 0, 10}, PlotRange -> All]

We can see that we get reasonably good list plot using the options MinRecursion and PrecisionGoal given to NIntegrate in step 1. (The grid lines show the sampling points.)

AbsoluteTiming[
 pres = PPartialIntegrations[Sin[#] + 2 Exp[-100 (# - 5)^2] &, {0, 10}, 
    MinRecursion -> 1, PrecisionGoal -> 10];
]
ListLinePlot[pres, ImageSize -> Large, PlotRange -> All, 
 GridLines -> {pres[[All, 1]], None}, 
 GridLinesStyle -> {Directive[GrayLevel[0.9]], None}]

(* {0.062035, Null} *)

Animation showing adaptivity

Following a misinterpreted suggestion by J.M. I made this animation that shows the effect the options of PPartialIntegrations.