74. Numerical integration#
Anne Kværnø, André Massing
Date: Jan 25, 2021
If you want to have a nicer theme for your jupyter notebook, download the cascade stylesheet file tma4125.css and execute the next cell:
from math import factorial
from IPython.core.display import HTML
from matplotlib.pyplot import *
from numpy import *
def css_styling():
try:
with open("tma4125.css") as f:
styles = f.read()
return HTML(styles)
except FileNotFoundError:
pass # Do nothing
# Comment out next line and execute this cell to restore the default notebook style
css_styling()
75. Introduction#
Imagine you want to compute the finite integral
The “usual” way is to find a primitive function
While there are many analytical integration techniques and extensive tables to determine definite integral for many integrands, more often than not it may not feasible or possible to compute the integral. For instance, what about
Finding the corresponding primitive is highly likely a hopeless
endeavor. And sometimes there even innocent looking functions
like
A numerical quadrature or a quadrature rule is a formula for
approximating such definite integrals
where
We sometimes might write
If the function
75.1. Example 1: Quadrature rules from previous math courses#
The trapezoidal rule, the midpoint rule and Simpson’s rule known from previous courses are all examples of numerical quadratures, and we quickly review them here:
Midpoint rule is the simplest possible quadrature rule defined by
The node is given by the midpoint,
Trapezoidal rule is given by
and thus the nodes are defined by
Finally, Simpson’s rule which you know from M1, is defined as follows:
which we identify as quadrature rule with 3 points
In this note we will see how quadrature rules can be constructed from integration of interpolation polynomials. We will demonstrate how to do error analysis and how to find error estimates. In the sequel, we will use material from Preliminaries, section 3.2, 4 and 5.
As usual, we import the necessary modules before we start computing.
%matplotlib inline
newparams = {
"figure.figsize": (8.0, 4.0),
"axes.grid": True,
"lines.markersize": 8,
"lines.linewidth": 2,
"font.size": 14,
}
rcParams.update(newparams)
76. Quadrature based on polynomial interpolation.#
This section relies on the content of the note on polynomial interpolation, in particular the section on Lagrange polynomials.
Choose
We will then use
with the cardinal functions
the following quadrature formula is obtained
where the weights in the quadrature are simply the integral of the cardinal functions over the interval
Let us derive three schemes for integration over the interval
76.1. Example 2: The trapezoidal rule revisited#
Let
and the corresponding quadrature rule is the trapezoidal rule (usually
denoted by
76.2. Example 3: Gauß-Legendre quadrature for #
Let $x_0=1/2 + \sqrt{3}/6$ and $x_1 = 1/2 - \sqrt{3}/6$. Then
The quadrature rule is
76.3. Example 4: Simpson’s rule revisited#
We construct Simpson's rule on the interval $[0,1]$ by choosing the nodes $x_0=0$, $x_1=0.5$ and $x_2=1$. The corresponding cardinal functions arewhich gives the weights
such that
76.4. Exercise 1: Accuracy of some quadrature rules#
Use the QR
function
def QR(f, xq, wq):
"""Computes an approximation of the integral f
for a given quadrature rule.
Input:
f: integrand
xq: quadrature nodes
wq: quadrature weights
"""
n = len(xq)
if n != len(wq):
raise RuntimeError("Error: Need same number of quadrature nodes and weights!")
return np.array(wq) @ f(np.array(xq))
to compute an approximate value of integral
using the quadrature rules from Example Example 2: The trapezoidal rule revisited-Example 4: Simpson’s rule revisited.
Tabulate the corresponding quadrature errors
# Insert your code here.
Hint.
You can start with (fill in values for any
# Define function
def f(x):
return ...
# Exact integral
int_f = ...
# Trapezoidal rule
xq = [...]
wq = [..., ...]
qr_f = QR(f, xq, wq)
print(f"Q[f] = {qr_f}")
print(f"I[f] - Q[f] = {int_f - qr_f:.10e}")
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[5], line 15
12 xq = [...]
13 wq = [..., ...]
---> 15 qr_f = QR(f, xq, wq)
16 print(f"Q[f] = {qr_f}")
17 print(f"I[f] - Q[f] = {int_f - qr_f:.10e}")
Cell In[3], line 12, in QR(f, xq, wq)
10 n = len(xq)
11 if n != len(wq):
---> 12 raise RuntimeError("Error: Need same number of quadrature nodes and weights!")
13 return np.array(wq) @ f(np.array(xq))
RuntimeError: Error: Need same number of quadrature nodes and weights!
76.4.1. Remarks#
We observe that with the same number of quadrature points, the Gauß-Legendre quadrature gives a much more accurate answer then the trapezoidal rule. So the choice of nodes clearly matters. Simpon’s rule gives very similar results to Gauß-Legendre quadrature, but it uses 3 instead of 2 quadrature nodes.
77. Degree of exactness and an estimate of the quadrature error#
Motivated by the previous examples, we know take a closer look at how assess the quality of a method. We start with the following definition.
77.1. Definition 1: The degree of exactness#
A numerical quadrature has degree of exactness
Since both integrals and quadratures are linear in the integrand
Observation:
All quadratures constructed from Lagrange interpolation polynomials in
77.2. Exercise 2: Degree of exactness for some quadrature rules#
Show that the trapezoidal and midpoint rule from Example Example 1: Quadrature rules from previous math courses is of precision 1
Show that the Gauß-Legendre quadrature for 2 points from Example Example 3: Gauß-Legendre quadrature for
is of precision 3.Show that Simpson’s rule from Example Example 4: Simpson’s rule revisited is also of precision 3.
# Insert your code here.
Hint. You can do this either using pen and paper (boring!) or numerically (more fun!), using the code from Example Exercise 1: Accuracy of some quadrature rules
Solution.
77.3. Error estimates#
77.4. Theorem 1: Error estimate for quadrature rule with degree of exactness #
Assume that
where
Proof.
Let
for some
which concludes the proof.
The previous theorem often The advantage of the previous theorem is that it is easy to prove. On downside is that the provided estimate can be rather crude, and often sharper estimates can be established. We give two examples here of some sharper estimates (but without proof).
77.5. Theorem 2: Error estimate for the trapezoidal rule#
For the trapezoidal rule, there is a
77.6. Theorem 3: Error estimate for Simpson’s rule#
For Simpson’s rule, there is a