Approximating optimal step size for numerical integrator

Question

I am trying to find a formula to approximate the ideal step size for the Trapezoid and Simpson's rules. As an example, consider the finite difference formula $$g(x,h) = \frac{f(x+h)-f(x-h)}{2h}$$ where $\lim_{h \rightarrow 0} g(x,h) = f'(x)$. The truncation error is bounded by $$\epsilon_{tr} \leq \frac{h^2}{6} f'''(x) $$ When computing the value of $g$ numerically, the value we get carries some rounding error, $$\tilde{g}(x,h) = \frac{f(x+h)(1\pm \epsilon_p) - f(x-h)(1\pm \epsilon_p)}{2h}$$ where $\epsilon_p$ is the "machine precision." The total rounding error is $$\epsilon_{ro} = \left|\tilde{g}(x,h) - g(x,h)\right| = \frac{\epsilon_p}{2h}\Bigl| \pm f(x+h) \pm f(x-h)\Bigr| \leq \frac{\epsilon_p}{h}\mathrm{max} \left|f(x)\right|$$ where the max is taken in the range $[x-h,x+h]$. If $h$ is very small we can approximate $$\epsilon_{ro} \approx \frac{\epsilon_p}{h}\left|f(x)\right|$$

Truncation error decreases as $h$ decreases while roundoff error increases as $h$ decreases. Finding the approximate value at which the two are of the same order of magnitude, we have $$\frac{\epsilon_p}{h_c}\left|f(x)\right| = \frac{h_c^2}{6} \left|f'''(x)\right|$$ or $$h_c \sim \left| \frac{6 \epsilon_p f(x)}{f'''(x)} \right|^{1/3}$$

(note: I'm not sure about the factor 6, as I have seen this quoted as a factor 3. I think this might be important, as we have two roundoff errors and they are as likely to combine as they are to cancel. See for example here)

Now I want to do something similar for the Trapezoidal and Simpson's rules. For the trapezoidal rule I have $$I_t(N) = \frac{h}{2}\left[f_0 + 2\sum_{j=1}^{N-1} f_j + f_{N}\right]$$ where the integration range is $[a,b]$, $f_n = f(x_n)$ and $x_n = a + (b-a)n/N$, and $h = (b-a)/N$. The truncation error is given by $$\epsilon_{tr} \leq \frac{h^2}{12}\left|f'(b) - f'(a)\right|$$ I find the roundoff error to be $$\epsilon_{ro} =\left|I_t(N) - \tilde{I}_t(N)\right| \sim \epsilon_p h \left|\sum_{j=0}^N (\pm f_j)\right|$$ Now I'm not really sure what to do, I could say that $$\epsilon_{ro} \leq \epsilon_p h \left|\sum_{j=0}^N \left|f_j\right| \right| < N \epsilon_p h \,\mathrm{max}\left|f(x)\right|$$ but this seems like a gross overestimation, in light of the point about the errors cancelling each other out. This would lead me to $$N_c \sim \sqrt{\frac{(b-a)}{12 \epsilon_p} \frac{f'(b) - f'(a)}{\mathrm{max}\left|f(x)\right|}}$$

Any help appreciated. One thing I am thinking is estimating the roundoff error as $$\epsilon_{ro} \sim \sqrt{N} \,\epsilon_p h\, \mathrm{avg}\left|f(x)\right|$$ where the $\sqrt{N}$ comes from thinking of the roundings as being a random walk and the average of the function could be easily estimated with just a few tens of points. This would give me $$N_c \sim \left| \frac{(b-a)}{12\epsilon_p} \frac{f'(b)-f'(a)}{\mathrm{avg}\left|f(x)\right|}\right|^{2/3}$$

I'm wondering if it can be improved.

Answer 1

I think you are missing an essential concept: We need an optimal stepsize for numerical differentiation because the condition number of numerical differentiation is infinite. For numerical integration, the condition number depends on $f$, but (for example) if $f$ is strictly positive, the problem is well-conditioned and the condition number of the quadrature sum is 1. You can basically compute an infinite amount of time with the trapezoidal rule and the approximation just gets better until $N\epsilon \approx \left\|f\right\|_1$, and even then you can use Kahan summation and compute until the end of the lifetime of the universe, the accuracy just getting better and better. You can learn the proof of these facts by reading Higham's classic paper The Accuracy of Floating Point Summation.