1. You are not fitting your data using sgolayfilt. All that does is try to reduce the noise in the data, assuming the underlying relationship is a smooth one, that is well represented locally by a low order polynomial.
2. Of course cumtrapz and trapz produce different results, since cumtrapz is a cumulative integral. So cumtrapz produces a vector result, whereas trapz will produce a scalar value, i/e/. the total area. HOWEVER, the final element of the vector returned by cumtrapz will be the same as that returned from trapz. The only reason that will not be true is if there is some amount of trash in the least significant bits of the result. That is something you cannot control. But there will be no significant difference otherwise.
Why do you need to crosscheck the result of trapz? This is code that has existed for dozens of years, written by professionals, been used millions, even billions of times. As much, it is trivial code to write. It is correct. I'm not worried. So what are you trying to verify? If you are still really worried, just test trapz against a known result.
Again, an sgolayfilt curve is not a curve fit by some model. There is no explicit curve that has been produced. The result is just a list of values at discrete points. So trapz is the reasonable solution. In fact, trapz in the presence of noise, trapz will be the lowest variance estimator of the integral. Other methods, for example, Simpson's rule, will have a slightly higher variance, although you may find a reduction in the bias of that estimate. Higher order Newton-Cotes rules will show progressively higher variances due to the noise still present in the signal, so at some point you will not gain by use of a increasingly higher order rule when applied to a noisy signal.
You might even be able to argue that the Newton-Cotes rule should NOT be a higher order rule than the implicit local polynomial order from the Savitsky-Golay filter. (That is something I'd need to think about.)
3. Could you use a Savitsky-Golay class filter to produce a direct estimate of the cumulative integral of the curve, rather than the function value at each point? It would seem so. Not even that difficult to write of course, but I don't know of any implementations on the file exchange.
