speed up integrating the same function over many overlapping intervals?

I need to compute many different integrals of the same function fn, essentially like this:

Xint = size(X);
for i=1:numel(X)
    Xint(i) = integral(fn,0,X(i));  % fn is continuous
end

The problem is that X is large and computation of fn is slow, so this simple loop takes quite a long time. My question is, what is the best way to speed up this loop?

It does seem like considerable speed-up should be possible, because each time through the loop the integral computation evaluates fn at many almost-the-same points from 0 to X(i)--the range of the X(i)'s is about 400-2000. But what's the best approach? E.g. maybe just order the X's and compute the integrals in pieces, 0 to min(X), then min(X) to 2nd-smallest-X, etc (but isn't there a danger that the errors of approximation will accumulate)? Or maybe precompute fn at a bunch of X's and use trapz (but how to choose those to-be-precomputed X's to ensure the desired level of accuracy)? Or precompute fn, spline it, and then let integral work with the spline_of_fn?

Given all the work that has been done on numerical methods of integration, it seems like someone must have studied how to do this most effectively, but I can't find anything on it. So, I'd appreciate any pointers or tips.

Thanks,

[Paul]

Might the best appoach depend on the properties of fn?

For this approach: "order the X's and compute the integrals in pieces," what is the concern of accumulation of errors? Are you concerned that

integral(fn,0,X(1)) + integral(fn,X(1),X(2));

has more error than

integral(fn,0,X(2)); %?

[Jeff Miller]

@Paul Yes, that is my concern, and even moreso as the sum accumulates more terms, e.g., out to

sum_at_X_end_minus_1 + integral(fn,X(end-1),X(end));

I am seeing big discrepancies between the sum accumulated piecewise like that versus

integral(fn,0,max(X))  %...and I assume this one is more accurate

[Paul]

Interesting. Do you have good reason to assume that the latter is more accurate? I wonder if this difference between the two approaches is unique to the function in question.

[Jeff Miller]

@Paul My reason for the assumption is not very good--it's just based on a vague understanding of how these integrals get computed to a convergence criterion, leaving a little residual error in each integral. Often I'd expect the piecewise errors to more or less cancel each other out when the pieces are summed, but I doubt that would be guaranteed for all functions.

It would be illuminating to experiment with full-range versus piece-wise integrals for some functions where the correct value can be computed analytically.

Answer 1

[Alan Stevens]

How about using cumtrapz to do the cumulative integral from 0 to max(X) just once, then select the sections specified by X(i).

[Jeff Miller]

Thanks for the suggestion. That is indeed much much faster. Thus far I have not been able to get satisfactory levels of accuracy, but that should presumably be possible by evaluating the function at a dense enough spacing of X values. (to be continued...)

[Walter Roberson]

integral() and the older functions such as quadgk() use adaptive step-sizes, with different methods for judging when adaption is needed, using different weights and number of evaluations, and using different methods for dealing with infinities.

"Black-box" functions are difficult to integrate reliably.

[Jeff Miller]

Thanks for your thoughts. Just a small followup comment:

The adaptive step size methods are no doubt very helpful in computing a single integral (say 0 to X1), but it seems all the computations involved in step size adaptation have to be redone from scratch when the integration range changes (say 0 to X2). My X vector can easily contain several hundred values, so computation of all the different integrals could in principle be speeded up by somehow saving some or all of the computations relevant to step size adaptation.

For comparison, I approximated a single integral using cumtrapz with a naive equal spacing of x values. Roughly, I needed so many x values (to match integral()'s precision) that the cumtrapz method took 2-3 times as long as integral(). But of course the cumtrapz method wins when the number of integrations increases, because integral() has to redo all the function calculations and cumtrapz 'remembers' them from the earlier computations. For repeated integrations, it seems like it should be possible in principle to have the best of both worlds--adaptive choice of step sizes but also some retention of the function values computed along the way. Not that I would have any idea how to code such a thing...

Answer 2

[Paul]

In your question you said that fn is continuous. If it's sufficiently smooth, ideally with continuous first derivative, you can use an ode solver to get the value of the integral at the points you specify in the tspan vector. For example:

fn = @(x) (sin(x));
odefun = @(t,x) (fn(t));
X = linspace(0,1,10)*2*pi;
[~,intvalues] = ode23(odefun,X,0); % note, X(1) = 0, using ode23 jsut to illustrate
plot(X,intvalues,X,-cos(X)+1,'o'),grid  % compare to truth solution

Of course you'll have to consider all the other options for the ode solver.

I have absolutely no idea how accurate or fast this solution is compared to any other solutions discussed in this thread, but at least it will give you the values of the integral at the specified values in X in one go.