Ch5.3: Gauss Quadrature

• Gauss quadrature achieves high accuracy for very few nodes.
  • Nodes and weights strategically computed.
  • Linear and quadratic interpolation (\( n=2 \), \( n=3 \))
• Greater accuracy and fewer flops than for Newton-Cotes.
  • Left figure: 2-Point Trapezoid Rule
  • Right figure: 2-Point Gauss Quadrature Rule

What do you call a farm that makes bad jokes?

Corny!

• Greek method for finding area of nonstandard region was to geometrically construct a square having same area.
  
• This is where the term quadrature comes from.

• On [a, b], a typical numerical integration method looks like

\[ \small{\int_{a}^b f(x)dx \cong \sum_{k=1}^N c_k f(x_k), \,\, c_k = \Delta x = \frac{b-a}{N} } \]

On [-1, 1] with two nodes, GQ uses non-endpoint trapezoid:

\[ \small{ \begin{aligned} c_1 &= c_2 = 1, \, \, x_1 = -\frac{1}{\sqrt{3}}, \, \, x_2 = \frac{1}{\sqrt{3}} \\ \int_{a}^b f(x)dx &\cong \sum_{k=1}^n c_k f(x_k) = 1 \cdot f(x_1) + 1 \cdot f(x_2) \end{aligned} } \]

For this example, \( f(x) = 1.7x^3 + 3x^2 + 2.4x +4 \)

www.desmos.com/calculator/udbnhupofh

For this example, \( f(x) = 1.7x^3 + 3x^2 + 2.4x +4 \)

www.desmos.com/calculator/uiuvdb8nsg

For \( n=2 \) data points, want GQ to be exact for polynomials of degree \( 2n-1 = 2(2)-1 = 3 \) or less on [-1,1]:

\[ \small{ \int_{-1}^1 P_3(x)dx = c_1 \cdot P_3(x_1) + 1 \cdot P_3(x_2) } \]

• For \( n=2 \) data points, Gauss quadrature is exact for polynomials of degree \( 2n-1 = 3 \) or less on [-1,1].
• For this example, \( f(x) = 1.7x^3 + 3x^2 + 2.4x +4 \).

• For \( N \) data points, Gauss quadrature is exact for polynomials of degree \( 2N-1 \) or less on [-1,1].
• The \( x_k \) are also the roots of Legendre polynomials.

gaussLegendre(2,0,1)

$x
[1] 0.2113249 0.7886751

$w
[1] 0.5 0.5

gaussLegendre(3,1,10)

$x
[1] 2.014315 5.500000 8.985685

$w
[1] 2.5 4.0 2.5

• In STAT 200, probabilities are calculated from integrals of standard normal distribution:

\[ \int_a^b e^{-x^2} dx \]

www.desmos.com/calculator/chjetmv2ae

gaussint <- function(f,x,w) {
  y <- f(x)
  return (sum(y*w))}

w <- c(1,1)
x <- c(-1/sqrt(3),1/sqrt(3))
f <- function(x){x^2}
gaussint(f,x,w)

[1] 0.6666667

\[ \small{ \int_{a}^b f(x)dx \cong c_1 \cdot f(x_1) + c_2 \cdot f(x_2)} \]

• To switch from [-1,1] to \( [a,b] \), need to use a \( u \)-substitution:

\[ \small{ u = \frac{2x-a-b}{b-a}, \,\, a \leq x \leq b} \]

• R program:

gaussquad2 <- function(f,a,b) {
  w <- c(1,1)
  x <- c(-1/sqrt(3),1/sqrt(3))
  d <- (b-a)/2
  s <- (b+a)/2
  u <- d*x + s
  y <- f(d*x + s)
  return(d*sum(y*w))}

gaussquad3 <- function(f,a,b) {
  w <- c(5/9,8/9,5/9)
  x <- c(-sqrt(3/5),0,sqrt(3/5))
  d <- (b-a)/2
  s <- (b+a)/2
  u <- d*x + s
  y <- f(d*x + s)
  return(d*sum(y*w))}