0
$\begingroup$

Consider an integral to evaluate:

$$I=\int_a^b f(x)\,\mathrm dx.\tag1$$

In the Monte-Carlo quadrature, if $I_n$ is $n$th estimate of $I$, then for $(n+1)$th estimate we have:

$$I_{n+1}=\frac{nI_n+f(\operatorname{rand}())(b-a)}{n+1}.\tag2$$

This lets one track current estimate at each step, with constant rate of updating.

Monte-Carlo quadrature converges as $\delta I_n=O(n^{-1/2})$, which is quite slow. A faster approach could be the trapezoidal rule. With it, if we define (for $n$ being number of points)

$$S_{n+1,\,\alpha,\,\beta}=\sum_{k=1}^{n-1+\beta} f\left(a+\frac{b-a}n (k+\alpha)\right),$$

then

$$\begin{align} I_{n+1} &=\frac{b-a}{n }\left(\frac{f(a)+f(b)}2+S_{n+1,\,0,\,0}\right),\\ I_{2n+1}&=\frac{I_{n+1}}2+\frac{b-a}{2n}S_{n+1,\,0.5,\,1}. \end{align} \tag3$$

Its convergence in the number of points taken is better: $\delta I_n=O(n^2)$. But to update without needless re-evaluation of $f$ and without storing most of its values (except $f(a)+f(b)$), we have to do twofold increase of $n$ each time. This means that rate of updating drops exponentially each update, although the precision does increase quadratically with number of points.

I'm looking for a quadrature method, which would combine the good properties of the two methods described above:

  1. Rapid convergence (quadratic in number of points or better),
  2. Reuse of previous result without storing most of the values of integrand,
  3. Possibility of constant-rate updates.

Is there any such method?

$\endgroup$
  • $\begingroup$ Well you could always do trapezoidal with $2$ steps, then $4$ steps, then $8$ etc. This allows you to use the previous computed values and you still have quadratic convergence. Wouldn't this satisfy all your requirements? $\endgroup$ – Winther Apr 7 at 19:58
  • $\begingroup$ @Winther this will not satisfy requirement #3: doing more steps means taking more time for an update. $\endgroup$ – Ruslan Apr 7 at 20:31
0
$\begingroup$

You could do a recursive method where on each interval you estimate the error (such as for each interval, do another function eval at the midpoint and compare trapezoid to simpson's). If the error is small enough there, then you can add the integral on the interval to the total integral and ignore it. If not, then subdivide the integral in two, and repeat.

You're still having to store some function values and there's a lot more bookkeeping, but the advantage of this method vs just halving the step size is that you won't have to do all the work in places where the function is nice.

That said, what I'm wondering is how natural properties 2 and 3 are, and if the spectral convergence of something like Gauss-Kronrod (https://en.wikipedia.org/wiki/Gauss–Kronrod_quadrature_formula) would be enough to change your mind.

$\endgroup$
  • $\begingroup$ I'm aware of Gauss-Kronrod, it's not suitable for my needs, where updating without increasing waiting times is crucial. It's even more important than rate of convergence, since currently I'm using Monte-Carlo instead of trapezoid rule. I just wanted to speed it up, without losing the important feature. $\endgroup$ – Ruslan Apr 8 at 14:21
  • $\begingroup$ Could you go into a little more detail about your use-case and the constraints thereof? $\endgroup$ – JCK Apr 8 at 14:24
  • $\begingroup$ It's a light scattering simulation, intended to make current progress visible by the user. Each frame a new integration point is added, and the user can see what the scene looks like. Then, once the scene looks good enough, the simulation can be stopped. But if the updates take longer and longer, the user soon won't have any patience to wait until next frame appears. $\endgroup$ – Ruslan Apr 8 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.