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

Why are there gates around cemeteries?

 Because people are dying to get in!

• 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)} \]

We can map between [-1,1] and \( [a,b] \) with a change of variable:

\[ \small{ \begin{aligned} u &= \frac{2x-a-b}{b-a}, \,\, a \leq x \leq b \\ x & = \frac{b-a}{2}u + \frac{b+a}{2}, \,\, -1 \leq u \leq 1 \ \end{aligned}} \]

For purposes of integration, we have

\[ \small{ \begin{aligned} x & = \frac{b-a}{2}u + \frac{b+a}{2}, \,\, -1 \leq u \leq 1 \\ dx &= \frac{b-a}{2}du \end{aligned}} \]

Then

\[ \small{ \begin{aligned} \int_a^b f(x)dx &= \frac{b-a}{2} \int_{-1}^1 f\left(\frac{b-a}{2} u+\frac{b+a}{2} \right)du \\ &\cong \frac{b-a}{2} \sum_{k=1}^n c_k \cdot f\left(\frac{b-a}{2} u_k+\frac{b+a}{2} \right) \end{aligned}} \]

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))}

\[ \small{ \begin{aligned} \int_a^b f(x)dx & \cong \frac{b-a}{2} \sum_{k=1}^n c_k f\left(\frac{b-a}{2} u_k+\frac{b+a}{2} \right) \end{aligned}} \]

\[ \small{ \begin{aligned} \int_0^1 x^2 dx &= \frac{1-0}{2} \int_{-1}^1 f\left(\frac{1-0}{2} u+\frac{1+0}{2} \right)du \\ & \cong \frac{1-0}{2} \sum_{k=1}^2 c_k \cdot \left(\frac{1-0}{2} u_k+\frac{1+0}{2} \right)^2 \\ &= \frac{1}{2} \cdot 1 \cdot \left(-\frac{1}{2\sqrt{3}} + \frac{1}{2} \right)^2 + \frac{1}{2} \cdot 1 \cdot \left(\frac{1}{2\sqrt{3}} +\frac{1}{2} \right)^2 \\ &= \frac{1}{8} \left[ \left(\frac{1}{3}-\frac{1}{\sqrt{3}} + 1\right) + \left( \frac{1}{3}+\frac{1}{\sqrt{3}} +1 \right)\right] = \frac{1}{8} \left(\frac{2}{3} + 2 \right) = \frac{1}{3} \end{aligned}} \]

Here we have used

\[ \small{ x = \frac{b-a}{2}u + \frac{b+a}{2}, \,\, dx = \frac{b-a}{2}du } \]

\[ \small{ \begin{aligned} \int_1^{1.5} e^{-x^2} dx &= \frac{1.5-1}{2} \int_{-1}^1 f\left(\frac{1.5-1}{2} u+\frac{1.5+1}{2} \right)du \\ & \cong \frac{1}{4} \sum_{k=1}^2 c_k \cdot e^{-\left(\frac{1}{4} u_k+\frac{5}{4} \right)^2} \\ &= \frac{1}{4} \cdot 1 \cdot e^{-\left(-\frac{1}{4\sqrt{3}} + \frac{5}{4} \right)^2} + \frac{1}{2} \cdot 1 \cdot e^{-\left(\frac{1}{4\sqrt{3}} +\frac{5}{4} \right)^2} \end{aligned}} \]

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))}