Today we talked about numerical integration. We covered these methods:
Both of the last two involve interpolation on each subinterval (with two nodes or three nodes, respectively).
Today we talked about the composite Simpson’s rule. For \(f \in C^4[a,b]\), the formula is
\[\int_a^b f(x)\,dx \approx \frac{h}{3}(f(x_0) + 4f(x_1) + 2f(x_2) + 4 f(x_3) + 2 f(x_4) + \ldots + 4f(x_{n-1}) + f(x_n)),\]
where \(h = \frac{b-a}{n}\) and \(x_k = a + k h\). This can be coded in Python as:
def simpson(f, a, b, n):
h = (b-a)/n
total = f(a)+f(b)
total += sum([4*f(a+k*h) for k in range(1,n,2)])
total += sum([2*f(a+k*h) for k in range(2,n,2)])
return h*total/3The error using the composite Simpson’s rule is \[\text{Error} = \frac{(a-b)}{180} h^4 f^{(4)}(\xi) ~ \text{ for some } ~ \xi \in [a,b].\]
We did the following examples in class.
Estimate \(\displaystyle\int_1^{10} \frac{1}{\sqrt{x^3+1}} \, dx\) using \(n = 100\).
Estimate \(\displaystyle\int_0^{\pi} \sin x \, dx\) using \(n = 100\). Then calculate the error (which is 2 minus the estimate, since \(\int_0^{\pi} \sin x \, dx = 2\)). What does the error formula say the worst case error should be?
How large should \(n\) be in order to estimate \(\displaystyle\int_1^2 \frac{1}{x} \, dx\) with an error less than \(10^{-12}\)?
One problem with the composite Simpson’s method is it can be hard to know how big to pick \(n\) to get a desired error. We did this example:
We discussed two ways to solve this problem. One way is to integrate the Taylor series for \(f\). That actually works pretty well here, but not every function has a nice Taylor series. The other way is to use Simpson’s rule, but then it is hard to estimate the error.
The error formula requires computing the fourth derivative, but this function doesn’t have a nice 4th derivative. There is a trick to avoid that computation.
Idea. If you double \(n\), then the error should get 16 times smaller (approximately). So if \(S_n\) denotes the composite Simpson’s method with \(n\) subintervals and \(E_n\) is the error term, then \[E_{n} \approx 16 E_{2n}\] If we calculate both \[S_n = \text{True Integral} + 16E_{2n}\] and \[S_{2n} = \text{True Integral} + E_{2n},\] then the difference should be about 15 times the error \(E_{2n}\). So \[E_{2n} \approx \frac{S_n - S_{2n}}{15}.\]
After that, we derived the following version of Stirling’s Formula by applying the composite trapezoid rule to the intergral \(\displaystyle\int_1^n \ln x \, dx\).
\[n! \sim e \sqrt{n} \left( \frac{n}{e} \right)^n.\]
You can find the details on the Wikipedia entry for Stirling’s formula.