Gaussian Quadrature with Legendre polynomials

I Legendre polynomials: P0 (x) = 1, P1 (x) = x,

(n + 1)Pn+1 (x) = (2n + 1)xPn (x) − nPn−1 (x) for n ≥ 1.
Gaussian Quadrature with Legendre polynomials
I Legendre polynomials: P0 (x) = 1, P1 (x) = x,
(n + 1)Pn+1 (x) = (2n + 1)xPn (x) − nPn−1 (x) for n ≥ 1.

I Gaussian quadrature
Z 1
f (x) dx ≈ c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn ),
−1

where x1 , x2 , · · · , xn ∈ (−1, 1) are distinct roots of Pn (x).

Gaussian Quadrature with Legendre polynomials
I Legendre polynomials: P0 (x) = 1, P1 (x) = x,
(n + 1)Pn+1 (x) = (2n + 1)xPn (x) − nPn−1 (x) for n ≥ 1.

I Gaussian quadrature
Z 1
f (x) dx ≈ c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn ),
−1

where x1 , x2 , · · · , xn ∈ (−1, 1) are distinct roots of Pn (x).

I DoP = 2n − 1 with choice

 
Z 1 Z 1 Y
def x − x j 
ci == Li (x) dx =  dx for i = 1, · · · , n.
−1 −1 xi − x j
j6=i
Gaussian Quadrature with Legendre polynomials
I Legendre polynomials: P0 (x) = 1, P1 (x) = x,
(n + 1)Pn+1 (x) = (2n + 1)xPn (x) − nPn−1 (x) for n ≥ 1.

I Gaussian quadrature
Z 1
f (x) dx ≈ c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn ),
−1

where x1 , x2 , · · · , xn ∈ (−1, 1) are distinct roots of Pn (x).

I DoP = 2n − 1 with choice

 
Z 1 Z 1 Y
def x − x j 
ci == Li (x) dx =  dx for i = 1, · · · , n.
−1 −1 xi − x j
j6=i

I Next: Estimate error in quadrature

Hermite Interpolation, with Legendre polynomial Pn (x)
I Given Legendre roots x1 , x2 , · · · , xn with
(x1 , f (x1 ), f 0 (x1 )), (x2 , f (x2 ), f 0 (x2 )), · · · , (xn , f (xn ), f 0 (xn )),
Hermite Interpolation, with Legendre polynomial Pn (x)
I Given Legendre roots x1 , x2 , · · · , xn with
(x1 , f (x1 ), f 0 (x1 )), (x2 , f (x2 ), f 0 (x2 )), · · · , (xn , f (xn ), f 0 (xn )),

I Interpolating polynomial H(x) of degree ≤ 2n − 1 satisfies

H(x1 ) = f (x1 ), H 0 (x1 ) = f 0 (x1 ),
H(x2 ) = f (x2 ), H 0 (x2 ) = f 0 (x2 ),
.. ..
. .
H(xn ) = f (xn ), H 0 (xn ) = f 0 (xn ).
Hermite Interpolation, with Legendre polynomial Pn (x)
I Given Legendre roots x1 , x2 , · · · , xn with
(x1 , f (x1 ), f 0 (x1 )), (x2 , f (x2 ), f 0 (x2 )), · · · , (xn , f (xn ), f 0 (xn )),

I Interpolating polynomial H(x) of degree ≤ 2n − 1 satisfies

H(x1 ) = f (x1 ), H 0 (x1 ) = f 0 (x1 ),
H(x2 ) = f (x2 ), H 0 (x2 ) = f 0 (x2 ),
.. ..
. .
H(xn ) = f (xn ), H 0 (xn ) = f 0 (xn ).

I Theorem: For each x ∈ [a, b], a number ξ(x) between

x1 , x2 , · · · , xn (hence ∈ (a, b)) exists with
f (2n) (ξ(x))
f (x) = H(x)+R(x), R(x) = (x−x1 )2 (x−x2 )2 · · · (x−xn )2 .
(2n)!
Hermite Interpolation, with Gaussian Quadrature

Z 1 Z 1 Z 1
f (x) dx = H(x) dx + R(x) dx
−1 −1 −1
Hermite Interpolation, with Gaussian Quadrature

Z 1 Z 1 Z 1
f (x) dx = H(x) dx + R(x) dx
−1 −1 −1
Z 1
deg(H)≤2n−1
======== c1 P(x1 ) + c2 P(x2 ) + · · · + cn P(xn ) + R(x) dx
−1
Hermite Interpolation, with Gaussian Quadrature

Z 1 Z 1 Z 1
f (x) dx = H(x) dx + R(x) dx
−1 −1 −1
Z 1
deg(H)≤2n−1
======== c1 P(x1 ) + c2 P(x2 ) + · · · + cn P(xn ) + R(x) dx
−1
magic
===== c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn )
Z 1 (2n)
f (ξ(x))
+ (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
−1 (2n)!
Hermite Interpolation, with Gaussian Quadrature

Z 1 Z 1 Z 1
f (x) dx = H(x) dx + R(x) dx
−1 −1 −1
Z 1
deg(H)≤2n−1
======== c1 P(x1 ) + c2 P(x2 ) + · · · + cn P(xn ) + R(x) dx
−1
magic
===== c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn )
Z 1 (2n)
f (ξ(x))
+ (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
−1 (2n)!
= c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn )
f (2n) (ξ) 1
Z
+ (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
(2n)! −1
Hermite Interpolation, with Gaussian Quadrature

Z 1 Z 1 Z 1
f (x) dx = H(x) dx + R(x) dx
−1 −1 −1
Z 1
deg(H)≤2n−1
======== c1 P(x1 ) + c2 P(x2 ) + · · · + cn P(xn ) + R(x) dx
−1
magic
===== c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn )
Z 1 (2n)
f (ξ(x))
+ (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
−1 (2n)!
= c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn )
f (2n) (ξ) 1
Z
+ (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
(2n)! −1
def
== c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn ) + R
Gaussian Quadrature Error Estimate

1
f (2n) (ξ)
Z
R = (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
(2n)! −1
Gaussian Quadrature Error Estimate

1
f (2n) (ξ)
Z
R = (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
(2n)! −1
!
22n 3
(n!) (n − 1)! 4−n f (2n) (ξ)
= f (2n) (ξ) = O .
(2n + 1)! (2n)! (2n − 1)! (2n)!
Gaussian Quadrature Error Estimate

1
f (2n) (ξ)
Z
R = (x − x1 )2 (x − x2 )2 · · · (x − xn )2 dx
(2n)! −1
!
22n 3
(n!) (n − 1)! 4−n f (2n) (ξ)
= f (2n) (ξ) = O .
(2n + 1)! (2n)! (2n − 1)! (2n)!

!
1 4−n f (2n) (ξ)
Z
f (x) dx = c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn ) + O .
−1 (2n)!

Rapid convergence for smooth functions

Composite Simpson’s Rule
b−a
(n = 2m, xj = a + j h, h = n ,0 ≤ j ≤ n)
 
b m−1 m
(b−a)h4 (4)
Z
h X X
f (x)dx = f (a) + 2 f (x2j ) + 4 f (x2j−1 ) + f (b)− f (µ)
a 3 180
j=1 j=1
Limitations of Gaussian Quadrature
Simpson/Trapezoidal:
I Composite rules:
I Adding more equi-spaced points.
I Romberg extrapolation:
I Obtaining higher order rules from lower order rules.
I Adaptive quadratures:
I Adding more points only when necessary.

Gaussian Quadrature:
I points different for different n.

Gaussian Quadrature good for given n,

not as good for given tolerance.
 RR
Double Integral Rf (x, y ) dA

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Double Integral = Integral of Integral

Z Z Z b Z d 
f (x, y ) dA = f (x, y ) dy dx
R a c
Z b
= g (x) dx,
a
Z d
def
where g (x) == f (x, y ) dy .
c

Approach

Rb
I Approximate a g (x) dx with quadrature.
I For any given xi , approximate g (xi ) with quadrature.
RR Rb Rd
R f (x, y ) dA = Rb a g (x) dx, g (x) = cf (x, y ) dy
I Approximate a g (x) dx with n-point quadrature:
Z b
g (x) dx = c1 g (x1 ) + c2 g (x2 ) + · · · + cn g (xn ) + R(g )
a
RR Rb Rd
R f (x, y ) dA = Rb a g (x) dx, g (x) = cf (x, y ) dy
I Approximate a g (x) dx with n-point quadrature:
Z b
g (x) dx = c1 g (x1 ) + c2 g (x2 ) + · · · + cn g (xn ) + R(g )
a

I For 1 ≤ i ≤ n, approximate g (xi ) with m-point quadrature:

Z d
c2 f (xi , y2 )+· · ·+b
f (xi , y ) dy = cb1 f (x1 , y1 )+b b (xi , ·)).
cm f (xi , ym )+R(f
c
RR Rb Rd
R f (x, y ) dA = Rb a g (x) dx, g (x) = cf (x, y ) dy
I Approximate a g (x) dx with n-point quadrature:
Z b
g (x) dx = c1 g (x1 ) + c2 g (x2 ) + · · · + cn g (xn ) + R(g )
a

I For 1 ≤ i ≤ n, approximate g (xi ) with m-point quadrature:

Z d
c2 f (xi , y2 )+· · ·+b
f (xi , y ) dy = cb1 f (x1 , y1 )+b b (xi , ·)).
cm f (xi , ym )+R(f
c

n
Z Z !
X
f (x, y ) dA = ci g (xi ) + R(g )
R i=1
   
Xn Xm
=  c i  b (xi , ·)) + R(g )
cbj f (xi , yj ) + R(f
i=1 j=1
RR Rb Rd
R f (x, y ) dA = Rb a g (x) dx, g (x) = cf (x, y ) dy
I Approximate a g (x) dx with n-point quadrature:
Z b
g (x) dx = c1 g (x1 ) + c2 g (x2 ) + · · · + cn g (xn ) + R(g )
a

I For 1 ≤ i ≤ n, approximate g (xi ) with m-point quadrature:

Z d
c2 f (xi , y2 )+· · ·+b
f (xi , y ) dy = cb1 f (x1 , y1 )+b b (xi , ·)).
cm f (xi , ym )+R(f
c

n
Z Z !
X
f (x, y ) dA = ci g (xi ) + R(g )
R i=1
   
Xn Xm
=  c i  b (xi , ·)) + R(g )
cbj f (xi , yj ) + R(f
i=1 j=1
 
n X
m n
!
X X
=  ci cbj f (xi , yj ) + b (xi , ·))
ci R(f + R(g ).
i=1 j=1 i=1
RR Rb Rd
R f (x, y ) dA = a g (x) dx, g (x) = cf (x, y ) dy

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Z b
g (x) dx = c1 g (x1 ) + c2 g (x2 ) + · · · + cn g (xn ) + R(g )
a
Z d
f (xi , y ) dy = cb1 f (x1 , y1 ) + cb2 f (xi , y2 ) + · · · + cbm f (xi , ym )
c
b (xi , ·))
+ R(f
RR Rb Rd
R f (x, y ) dA = a g (x) dx, g (x) = cf (x, y ) dy

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Z b
g (x) dx = c1 g (x1 ) + c2 g (x2 ) + · · · + cn g (xn ) + R(g )
a
Z d
f (xi , y ) dy = cb1 f (x1 , y1 ) + cb2 f (xi , y2 ) + · · · + cbm f (xi , ym )
c
b (xi , ·))
+ R(f
 
n X
m n
Z Z !
X X
f (x, y ) dA =  ci cbj f (xi , yj ) + b (xi , ·))
ci R(f + R(g )
R i=1 j=1 i=1
n X
X m
≈ ci cbj f (xi , yj ).
i=1 j=1

I Double integral quadrature is a double sum.

I Need to work out total error for any given quadrature.
RR Pn Pm Pn b (xi , ·)) + R(g )
R
f (x, y ) dA = i=1 j=1 ci c
bj f (xi , yj ) + i=1 ci R(f

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Example, Simpson’s Rule with m = n = 3:
(x1 , x2 , x3 ) = a, a+b


I Simpson’s Rule on [a, b]: 2 , b .
c+d


I Simpson’s Rule on [c, d]: (y1 , y2 , y3 ) = c, 2 , d .
RR Pn Pm Pn b (xi , ·)) + R(g )
R
f (x, y ) dA = i=1 j=1 ci c
bj f (xi , yj ) + i=1 ci R(f

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Example, Simpson’s Rule with m = n = 3:
(x1 , x2 , x3 ) = a, a+b


I Simpson’s Rule on [a, b]: 2 , b .
c+d


I Simpson’s Rule on [c, d]: (y1 , y2 , y3 ) = c, 2 , d .
RR Pn Pm Pn b (xi , ·)) + R(g )
R
f (x, y ) dA = i=1 j=1 ci c
bj f (xi , yj ) + i=1 ci R(f

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Example, m = n = 3:
(x1 , x2 , x3 ) = a, a+b


I Simpson’s Rule on [a, b]: 2 ,b .

b
h5 (4)
Z
g (x) dx = c1 g (x1 ) + c2 g (x2 ) + c3 g (x3 ) − g (ξ),
a 90

b−a h
h= , (c1 , c2 , c3 ) = (1, 4, 1) .
2 3
(y1 , y2 , y3 ) = c, c+d


I Simpson’s Rule on [c, d]: 2 ,d .

d
k 5 ∂4f
Z
f (xi , y ) dy = cb1 f (xi , y1 )+b
c2 f (xi , y2 )+b
c3 f (xi , y3 )− (xi , ηi ),
c 90 ∂ 4 y

d −c k
k= , (b
c1 , cb2 , cb3 ) = (1, 4, 1) .
2 3
  
n X
m m
!
k5 ∂4f
Z Z X X
f (x, y ) dA =  ci cbj f (xi , yj ) − ci 4 (xi , ηi )
R 90 ∂ y
i=1 j=1 i=1

h5 b ∂ 4 f
Z
− (ξ, y )dy
90 a ∂ 4 x
  
n X
m m
!
Z Z X k 5 X ∂ 4f
f (x, y ) dA =  ci cbj f (xi , yj ) − ci 4 (xi , ηi )
R 90 ∂ y
i=1 j=1 i=1

h5 b ∂ 4 f
Z
− (ξ, y )dy
90 a ∂ 4 x
 
n X m m
!
X k 5 X ∂4f b
=  ci cbj f (xi , yj ) − ci (ξ, ηb)
90 ∂4y
i=1 j=1 i=1

h5 ∂4f
− (b − a) 4 (ξ, η)
90 ∂ x
  
n X
m m
!
Z Z X k 5 X ∂ 4f
f (x, y ) dA =  ci cbj f (xi , yj ) − ci 4 (xi , ηi )
R 90 ∂ y
i=1 j=1 i=1

h5 b ∂ 4 f
Z
− (ξ, y )dy
90 a ∂ 4 x
 
n X m m
!
X k 5 X ∂4f b
=  ci cbj f (xi , yj ) − ci (ξ, ηb)
90 ∂4y
i=1 j=1 i=1

h5 ∂4f
− (b − a) 4 (ξ, η)
90 ∂ x
n X
X m
= ci cbj f (xi , yj )
i=1 j=1
4 4
 
(b − a)(d − c) 4∂ f b 4∂ f
− k 4 (ξ, ηb) + h 4 (ξ, η) .
180 ∂ y ∂ x
RR Pn Pm Pn b (xi , ·)) + R(g )
R
f (x, y ) dA = i=1 j=1 ci c
bj f (xi , yj ) + i=1 ci R(f

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Example: Composite Simpson Rules.
RR Pn Pm Pn b (xi , ·)) + R(g )
R
f (x, y ) dA = i=1 j=1 ci c
bj f (xi , yj ) + i=1 ci R(f

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Example: Composite Simpson Rules.
I Composite Simpson on [a, b], xi =a+(i−1)h, 1≤i ≤n, h= b−a .
n−1

(b − a)h4 d ∂ 4 f (b − a)(d − c)h4 ∂ 4 f

Z
R(g ) = − 4
(ξ, y ) dy = − (ξ, η).
180 c ∂ x 180 ∂4x
RR Pn Pm Pn b (xi , ·)) + R(g )
R
f (x, y ) dA = i=1 j=1 ci c
bj f (xi , yj ) + i=1 ci R(f

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Example: Composite Simpson Rules.
I Composite Simpson on [a, b], xi =a+(i−1)h, 1≤i ≤n, h= b−a .
n−1

(b − a)h4 d ∂ 4 f (b − a)(d − c)h4 ∂ 4 f

Z
R(g ) = − 4
(ξ, y ) dy = − (ξ, η).
180 c ∂ x 180 ∂4x

d−c
I Composite Simpson on [c, d], yj =c+(j−1)k, 1≤ j ≤m, k = m−1 .
4 4
b (xi , ·)) = − (d − c)k ∂ f (xi , ηi )
R(f
180 ∂4y
n n
!
X (d − c)k 4 X ∂4f b
b (xi , ·)) = −
ci R(f ci (ξ, ηb)
180 ∂4y
i=1 i=1
(b − a)(d − c)k 4 ∂ 4 f b
= − (ξ, ηb).
180 ∂4y
Error Estimate, Double Integral with Composite Simpson

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }

Z Z n X
X m
f (x, y ) dA = ci cbj f (xi , yj )
R i=1 j=1
4f 4
 
(b − a)(d − c) 4∂ 4∂ f
− k 4 (ξ, ηb) + h 4 (ξ, η) ,
b
180 ∂ y ∂ x

b−a d −c
h= , k= .
n−1 m−1
 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

R = {(x, y ) | 1.2 ≤ x ≤ 2.4, 0.2 ≤ y ≤ 1. }

 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

R = {(x, y ) | 1.2 ≤ x ≤ 2.4, 0.2 ≤ y ≤ 1, }

 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

R = {(x, y ) | 1.2 ≤ x ≤ 2.4, 0.2 ≤ y ≤ 1, }

2.4 − 1.2 1 − 0.2
h = = 0.2, k = = 0.2.
7−1 5−1
 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

R = {(x, y ) | 1.2 ≤ x ≤ 2.4, 0.2 ≤ y ≤ 1, }

2.4 − 1.2 1 − 0.2
h = = 0.2, k = = 0.2.
7−1 5−1
∂4f 6 ∂4f 96
4
= − 4
, 4
=− .
∂ x (x + 2y ) ∂ y (x + 2y )4
∂4f 6
4
≤ ≈ 0.91553 for (x, y ) ∈ R,
∂ x (1.2 + 2 × 0.2)4
∂4f 96
≤ ≈ 14.648.
∂4y (1.2 + 2 × 0.2)4

(b − a)(d − c) 4 ∂ 4 f b ∂4f
Quad Error = k 4 (ξ, ηb) + h4 4 (ξ, η)
180 ∂ y ∂ x
(2.4 − 1.2)(1 − 0.2)
0.24 × 14.648 + 0.24 × 0.91553


≤
180
 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

R = {(x, y ) | 1.2 ≤ x ≤ 2.4, 0.2 ≤ y ≤ 1, }

2.4 − 1.2 1 − 0.2
h = = 0.2, k = = 0.2.
7−1 5−1
∂4f 6 ∂4f 96
4
= − 4
, 4
=− .
∂ x (x + 2y ) ∂ y (x + 2y )4
∂4f 6
4
≤ ≈ 0.91553 for (x, y ) ∈ R,
∂ x (1.2 + 2 × 0.2)4
∂4f 96
≤ ≈ 14.648.
∂4y (1.2 + 2 × 0.2)4

(b − a)(d − c) 4 ∂ 4 f b ∂4f
Quad Error = k 4 (ξ, ηb) + h4 4 (ξ, η)
180 ∂ y ∂ x
(2.4 − 1.2)(1 − 0.2)
0.24 × 14.648 + 0.24 × 0.91553


≤
180
≈ 1.328 × 10−4 .
 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

Z Z
log(x + 2y ) dA = 1.0360481 · · ·
R
 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

Z Z
log(x + 2y ) dA = 1.0360481 · · ·
R
Z Z 7 X
X 5
log(x + 2y ) dA ≈ ci cbj f (xi , yj )
R i=1 j=1
= 1.0360327 · · · .
 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

Z Z
log(x + 2y ) dA = 1.0360481 · · ·
R
Z Z 7 X
X 5
log(x + 2y ) dA ≈ ci cbj f (xi , yj )
R i=1 j=1
= 1.0360327 · · · .

Therefore
 
Z Z X7 X
5
log(x + 2y ) dA −  ci cbj f (xi , yj ) = 1.546 × 10−5
R i=1 j=1
 RR
Example: R log(x + 2y ) dA with n = 7, m = 5

Z Z
log(x + 2y ) dA = 1.0360481 · · ·
R
Z Z 7 X
X 5
log(x + 2y ) dA ≈ ci cbj f (xi , yj )
R i=1 j=1
= 1.0360327 · · · .

Therefore
 
Z Z X7 X
5
log(x + 2y ) dA −  ci cbj f (xi , yj ) = 1.546 × 10−5
R i=1 j=1

< 1.328 × 10−4 .

(Error Estimate)
2 Dimensional Gaussian Quadratures

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Z Z Z b Z d 
f (x, y ) dA = f (x, y ) dy dx.
R a c
2 Dimensional Gaussian Quadratures

R = {(x, y ) | a ≤ x ≤ b, c ≤ y ≤ d. }
Z Z Z b Z d 
f (x, y ) dA = f (x, y ) dy dx.
R a c

Perform change of variables

a+b b−a c +d d −c
x= + u, y= + v for u, v ∈ [−1, 1].
2 2 2 2
Double integral becomes
1
(b − a)(d − c)
Z Z Z
f (x, y ) dA = gb(u)du, where
R 4 −1

1  
a+b b−a c +d d −c
Z
def
gb(u) == f + u, + v dv .
−1 2 2 2 2
RR (b−a)(d−c) R 1
R f (x, y ) dA = 4 −1 g (u)du
b
1  
a+b b−a c +d d −c
Z
gb(u) = f + u, + v dv .
−1 2 2 2 2
R1
I n-point Gaussian quadrature for −1 gb(u)du:
Z 1
gb(u)du ≈ c1 gb(u1 ) + c2 gb(u2 ) + · · · + cn gb(un ).
−1
RR (b−a)(d−c) R 1
R f (x, y ) dA = 4 −1 g (u)du
b
1  
a+b b−a c +d d −c
Z
gb(u) = f + u, + v dv .
−1 2 2 2 2
R1
I n-point Gaussian quadrature for −1 gb(u)du:
Z 1
gb(u)du ≈ c1 gb(u1 ) + c2 gb(u2 ) + · · · + cn gb(un ).
−1

I For 1 ≤ i ≤ n, let
1  
a+b b−a d −c
Z
c +d
xi = + ui , gb(ui ) = f xi , + v dv .
2 2 −1 2 2
RR (b−a)(d−c) R 1
R f (x, y ) dA = 4 −1 g (u)du
b
1  
a+b b−a c +d d −c
Z
gb(u) = f + u, + v dv .
−1 2 2 2 2
R1
I n-point Gaussian quadrature for −1 gb(u)du:
Z 1
gb(u)du ≈ c1 gb(u1 ) + c2 gb(u2 ) + · · · + cn gb(un ).
−1

I For 1 ≤ i ≤ n, let
1  
a+b b−a d −c
Z
c +d
xi = + ui , gb(ui ) = f xi , + v dv .
2 2 −1 2 2

I m-point Gaussian quadrature for gb(ui ),

   
c +d d −c c +d d −c
gb(ui ) ≈ cb1 f xi , + v1 + · · · + cbm f xi , + vm
2 2 2 2
def c + d d −c
= cb1 f (xi , y1 ) + · · · + cbm f (xi , ym ) , yj = + vj .
2 2
RR (b−a)(d−c) R 1
R f (x, y ) dA = 4 −1 g (u)du
b

Z 1
gb(u)du ≈ c1 gb(u1 ) + c2 gb(u2 ) + · · · + cn gb(un )
−1
gb(ui ) ≈ cb1 f (xi , y1 ) + · · · + cbm f (xi , ym ) .

So we have Gaussian quadrature for double integral:

n m
(b − a)(d − c) X X
Z Z
f (x, y ) dA ≈ ci cbj f (xi , yj ) .
R 4
i=1 j=1
 RR
Example: R log(x + 2y ) dA = 1.036 · · · , n = 7, m = 5
I Gaussian quadrature approximation
n m
(b − a)(d − c) X X
Z Z
f (x, y ) dA ≈ ci cbj f (xi , yj ) ≈ 1.03604817065 · · ·
R 4
i=1 j=1

n m
(b − a)(d − c) X X
Z Z
f (x, y ) dA − ci cbj f (xi , yj ) ≈ 6.4×10−10 .
R 4
i=1 j=1
 RR
Example: R log(x + 2y ) dA = 1.036 · · · , n = 7, m = 5
I Gaussian quadrature approximation
n m
(b − a)(d − c) X X
Z Z
f (x, y ) dA ≈ ci cbj f (xi , yj ) ≈ 1.03604817065 · · ·
R 4
i=1 j=1

n m
(b − a)(d − c) X X
Z Z
f (x, y ) dA − ci cbj f (xi , yj ) ≈ 6.4×10−10 .
R 4
i=1 j=1

I Simpson rule approximation

Z Z n X
X m
f (x, y ) dA ≈ ci cbj f (xi , yj ) ≈ 1.03603270963 · · · .
R i=1 j=1

Z Z n X
X m
f (x, y ) dA − ci cbj f (xi , yj ) ≈ 1.5 × 10−5 .
R i=1 j=1
 RR
Example: R log(x + 2y ) dA = 1.036 · · · , n = 7, m = 5
I Gaussian quadrature approximation
n m
(b − a)(d − c) X X
Z Z
f (x, y ) dA ≈ ci cbj f (xi , yj ) ≈ 1.03604817065 · · ·
R 4
i=1 j=1

n m
(b − a)(d − c) X X
Z Z
f (x, y ) dA − ci cbj f (xi , yj ) ≈ 6.4×10−10 .
R 4
i=1 j=1

I Simpson rule approximation

Z Z n X
X m
f (x, y ) dA ≈ ci cbj f (xi , yj ) ≈ 1.03603270963 · · · .
R i=1 j=1

Z Z n X
X m
f (x, y ) dA − ci cbj f (xi , yj ) ≈ 1.5 × 10−5 .
R i=1 j=1

Gaussian quadrature much more accurate .

Gaussian Quadrature: how good is it?

I Gaussian quadrature
Z 1
f (x) dx = c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn )
−1
!
4−n f (2n) (ξ)
+O
(2n)!
≈ c1 f (x1 ) + c2 f (x2 ) + · · · + cn f (xn ),

where x1 , x2 , · · · , xn ∈ (−1, 1) are distinct roots of Legendre

Polynomial Pn (x).
I Error tiny for large n and f (2n) (ξ) = O(1).

Improper Integral: what if maxx∈[−1,1] f (2n) (x)  1?

 R1
Gaussian Quadrature for √ 1
−1 f (x)dx, f (x) = 1−x 2

n Quadrature Value
10 2.9758
50 3.1071
250 3.1346
1000 3.1399
 R1
Gaussian Quadrature for √ 1
−1 f (x)dx, f (x) = 1−x 2

n Quadrature Value
10 2.9758
50 3.1071
250 3.1346
1000 3.1399

I f (x) = √ 1 not smooth in [−1.1]

1−x 2
 R1
Gaussian Quadrature for √ 1
−1 f (x)dx, f (x) = 1−x 2

n Quadrature Value
10 2.9758
50 3.1071
250 3.1346
1000 3.1399

I f (x) = √ 1 not smooth in [−1.1]

1−x 2

π
Z 1 Z
dx x=sinθ 2
√ ==== 1 dθ = 3.1416 · · ·
−1 1 − x2 − π2
 R1 f (x)
Improper Integral: −1
√
1−x 2
dx
I Change of variable:
π π
x = sin θ, θ ∈ [− , ].
2 2
I p p
dx = cos θ, 1 − x 2 = 1 − sin2 θ = cos θ.
I New integral becomes proper.
Z 1 Z π
f (x) 2
√ dx = f (sin θ) dθ.
1−x 2
−1 − π2
 R1 f (x)
Improper Integral: −1
√
1−x 2
dx
I Change of variable:
π π
x = sin θ, θ ∈ [− , ].
2 2
I p p
dx = cos θ, 1 − x 2 = 1 − sin2 θ = cos θ.
I New integral becomes proper.
Z 1 Z π
f (x) 2
√ dx = f (sin θ) dθ.
1−x 2
−1 − π2

I Example: for f (x) = x 2 , with 20-point Gaussian quadrature:

π
1
x2
Z Z
2
√ dx = sin2 θ dθ ≈ 1.57079632679490,
−1 1 − x2 − π2

accurate to 15 digits.
 Rb g (x)
Improper Integral: a f (x)dx, f (x) = (x−a)p , for p < 1

I Integral is improper if g (a) 6= 0:

g (x)
lim = ∞.
x→a (x − a)p
 Rb g (x)
Improper Integral: a f (x)dx, f (x) = (x−a)p , for p < 1

I Integral is improper if g (a) 6= 0:

g (x)
lim = ∞.
x→a (x − a)p

I Integral is not defined for p ≥ 1:

Z b Z b
1 1
p
dx = lim + p
dx
a (x − a) M→a M (x − a)
(x − a)1−p b (b − a)1−p
= lim + M = .
M→a 1−p 1−p
 Rb g (x)
Proper Method for improper Integral a (x−a)p dx
I Assume a Taylor expansion on g (x):
g 00 (a) g (n)
g (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (a)(x−a)n +· · · .
2 n!
 Rb g (x)
Proper Method for improper Integral a (x−a)p dx
I Assume a Taylor expansion on g (x):
g 00 (a) g (n)
g (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (a)(x−a)n +· · · .
2 n!

I Choose a (k + 1)-term approximation:

g 00 (a) g (k) (a)
Pk (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (x−a)k .
2 k!
 Rb g (x)
Proper Method for improper Integral a (x−a)p dx
I Assume a Taylor expansion on g (x):
g 00 (a) g (n)
g (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (a)(x−a)n +· · · .
2 n!

I Choose a (k + 1)-term approximation:

g 00 (a) g (k) (a)
Pk (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (x−a)k .
2 k!

I Improper integral decomposition

Z b Z b Z b
g (x) Pk (x) g (x) − Pk (x)
p
dx = p
dx + dx.
a (x − a) a (x − a) a (x − a)p
 Rb g (x)
Proper Method for improper Integral a (x−a)p dx
I Assume a Taylor expansion on g (x):
g 00 (a) g (n)
g (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (a)(x−a)n +· · · .
2 n!

I Choose a (k + 1)-term approximation:

g 00 (a) g (k) (a)
Pk (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (x−a)k .
2 k!

I Improper integral decomposition

Z b Z b Z b
g (x) Pk (x) g (x) − Pk (x)
p
dx = p
dx + dx.
a (x − a) a (x − a) a (x − a)p

I First integral is a simple sum:

b k Z b (j) k
g (a) (x − a)j g (j) (a)
Z
Pk (x) X X
dx = dx = (b−a)j+1−p .
a (x − a)p a j! (x − a) p j!(j + 1 − p)
j=0 j=0
 Rb g (x)
Proper Method for improper Integral a (x−a)p dx
I Assume a Taylor expansion on g (x):

g 00 (a) g (n)
g (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (a)(x−a)n +· · · .
2 n!
I Choose a (k + 1)-term approximation:

g 00 (a) g (k) (a)

Pk (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (x−a)k .
2 k!
I Improper integral decomposition
Z b Z b Z b
g (x) Pk (x) g (x) − Pk (x)
p
dx = p
dx + dx.
a (x − a) a (x − a) a (x − a)p
I Second integral is a proper integral for k  p:
 
b Z b ∞ (j) (a)
g (x) − Pk (x)
Z X g
dx = (x−a)k+1−p  (x − a)j−k−1  dx.
a (x − a)p a j!
j=k+1
Summary of Proper Method
I Choose a (k + 1)-term approximation:

g 00 (a) g (k) (a)

Pk (x) = g (a)+g 0 (a)(x−a)+ (x−a)2 +· · ·+ (x−a)k .
2 k!
I Improper integral decomposition
b Z b Z b
g (x) − Pk (x)
Z
g (x) Pk (x)
p
dx = p
dx + dx.
a (x − a) (x − a) (x − a)p
a a

k Z b
X g (j) (a) j+1−p 
=  (b − a) + G (x) dx,
j!(j + 1 − p) a
j=0
(
0, if x = a,
where G (x) = g (x)−Pk (x)
(x−a)p . if x > a.

Only require g (a), g 0 (a), · · · , g (k) (a)

 R1ex
Example: 0
√
x
dx
I Take a 5-term Taylor expansion

x2 x3 x4
P4 (x) = 1 + x + + + .
2 6 24
I
!
1 1 √ x 3/2 x 5/2 x 7/2
Z Z
P4 (x) 1
√ dx = √ + x+ + + dx
0 x 0 x 2 6 24
2 1 1 1
= 2+ + + +
3 5 21 108
≈ 2.9235450.

I Composite Simpson’s rule with n = 4, h = n1 = 0.25 on

Z 1 (
0, if x = 0,
G (x) dx, where G (x) = e x −P4 (x)
0
√
x
. if x > 0.
 R1 ex
Example: 0
√
x
dx
x G (x)
0.0 0
0.25 0.0000170
.
0.50 0.0004013
0.75 0.0026026
1.0 0.0099485

I Composite Simpson rule, with n = 4, h = n1 = 0.25

Z 1
0.25
G (x)dx ≈ (0 + 4 × 0.0000170 + 2 × 0.0004013
0 3
+ 4 × 0.0026026 + 0.0099485) ≈ 0.0017691.
 R1 ex
Example: 0
√
x
dx
x G (x)
0.0 0
0.25 0.0000170
.
0.50 0.0004013
0.75 0.0026026
1.0 0.0099485

I Composite Simpson rule, with n = 4, h = n1 = 0.25

Z 1
0.25
G (x)dx ≈ (0 + 4 × 0.0000170 + 2 × 0.0004013
0 3
+ 4 × 0.0026026 + 0.0099485) ≈ 0.0017691.

I Improper integral decomposition

Z 1 x Z 1 Z 1 x
e P4 (x) e − P4 (x)
√ dx = √ dx + √ dx
0 x 0 x 0 x
≈ 2.9235450 + 0.0017691 = 2.9253141.
 Rb g (x)
Improper Integral: a f (x)dx, f (x) = (b−x)p , for p < 1
I Change of variable: z = −x
I Left endpoint improper Integral
Z b Z −a
g (x) g (−z)
dx = dz
a (b − x)p −b (z − (−b))p
 R∞ g (x)
Improper Integral: a f (x)dx, f (x) = xp , for p > 1
I Integral is not defined for p ≤ 1:
Z ∞ Z M
1 1
dx = lim dx
a xp M→∞ a xp
x 1−p a1−p
= lim |∞
a = .
M→∞ 1−p p−1

I In general, change of variable z = x1 , assuming a > 0.

dz
dx = − .
z2
I Left end improper integral:
1
Z ∞ Z
g (x) a 1
dx = g ( ) z p−2 dz.
a xp 0 z
 R∞ sin( x1 )
Example: 1 f (x)dx, f (x) = x 3/2
I Change of variable z = x1 ,:
sin( x1 )
Z ∞ Z 1 Z 1
3/2 −2 sin(z)
3/2
dx = sin(z) z z dz = √ dz.
1 x 0 0 z
 R∞ sin( x1 )
Example: 1 f (x)dx, f (x) = x 3/2
I Change of variable z = x1 ,:
sin( x1 )
Z ∞ Z 1 Z 1
3/2 −2 sin(z)
3/2
dx = sin(z) z z dz = √ dz.
1 x 0 0 z

z3
I Choose 5-term Taylor expansion on sin(z): P4 (z) = z − 6.
 R∞ sin( x1 )
Example: 1 f (x)dx, f (x) = x 3/2
I Change of variable z = x1 ,:
sin( x1 )
Z ∞ Z 1 Z 1
3/2 −2 sin(z)
3/2
dx = sin(z) z z dz = √ dz.
1 x 0 0 z

z3
I Choose 5-term Taylor expansion on sin(z): P4 (z) = z − 6.
I
 
z3
Z 1
sin(z)
Z 1
z− z3 sin(z) −
Z z
1 − 6
√ dz 6
√ dz + √
= dz
0 z 0 z 0 z
 
Z 1 sin(z) − z − z 3
6
≈ 0.61904761 + √ dz.
0 z
I Composite Simpson’s rule with n = 16:
 
Z 1 sin(z) − z − z 3
6
√ dz ≈ 0.0014890097.
0 z
 R∞ sin( x1 )
Example: 1 f (x)dx, f (x) = x 3/2

I
 
z3
Z 1
sin(z)
Z
z−1 z3 sin(z) −Z z1− 6
√ dz 6
√ dz + √
= dz
0 z 0 z 0 z
 
Z 1 sin(z) − z − z 3
6
≈ 0.61904761 + √ dz
0 z
≈ 0.61904761 + 0.0014890097
= 0.62053661,

accurate to 8-digit.
Initial Value ODE
I The motion of a swinging pendulum

I Initial Value conditions

θ(t0 ) = θ0 , and θ0 (t0 ) = θ00 .

Initial Value ODE
I The motion of a swinging pendulum

I Initial Value conditions

θ(t0 ) = θ0 , and θ0 (t0 ) = θ00 .

When does ODE have a solution?

Initial Value ODE
I The motion of a swinging pendulum

I Initial Value conditions

θ(t0 ) = θ0 , and θ0 (t0 ) = θ00 .

When does ODE have a solution? How to compute it?

Lipschitz condition
Definition: function f (t, y ) satisfies a Lipschitz condition in the
variable y on a set D ⊂ R2 if a constant L > 0 exists with

|f (t, y1 ) − f (t, y2 )| ≤ L |y1 − y2 | ,

whenever (t, y1 ), (t, y2 ) are in D. L is Lipschitz constant.

I Example 1: Show that f (t, y ) = t|y | satisfies a Lipschitz
condition on the region

D = {(t, y ) | 0 ≤ t ≤ T }.

Solution: For any (t, y1 ), (t, y2 ) in D,

|f (t, y1 ) − f (t, y2 )| = |t|y1 | − t|y2 || ≤ t |y1 − y2 | ≤ L |y1 − y2 | ,

for L = T .
Lipschitz condition
Definition: function f (t, y ) satisfies a Lipschitz condition in the
variable y on a set D ⊂ R2 if a constant L > 0 exists with

|f (t, y1 ) − f (t, y2 )| ≤ L |y1 − y2 | ,

whenever (t, y1 ), (t, y2 ) are in D. L is Lipschitz constant.

I Example 2: Show that f (t, y ) = t y 2 does not satisfy any
Lipschitz condition on the region

D = {(t, y ) | 0 ≤ t ≤ T }.

Solution: Choose (T , y1 ), (T , y2 ) in D with y1 = 0, y2 > 0,

|f (T , y1 ) − f (T , y2 )|
= T y2 ,
|y1 − y2 |

which can be larger than L for any fixed L > 0.

Convex Set
Definition: A set D ⊂ R2 is convex if

whenever (t1 , y1 ) and (t2 , y2 ) ∈ D

−→ line segment (1 − λ) (t1 , y1 ) + +λ (t2 , y2 ) ∈ D for all λ ∈ [0, 1].
Theorem: Suppose f (t, y ) is defined on a convex set D ⊂ R2 . If a
constant L > 0 exists with
∂f
(t, y ) ≤ L, for all (t, y ) ∈ D,
∂y

then f satisfies a Lipschitz condition with Lipschitz constant L.

Theorem: Suppose f (t, y ) is defined on a convex set D ⊂ R2 . If a
constant L > 0 exists with
∂f
(t, y ) ≤ L, for all (t, y ) ∈ D,
∂y

then f satisfies a Lipschitz condition with Lipschitz constant L.

I Example 1: Show that f (t, y ) = t y 2 satisfies Lipschitz
condition on the region

D = {(t, y ) | 0 ≤ t ≤ T, −Y ≤ y ≤ Y } .

Solution:
∂f ∂f
(t, y ) = 2ty , (t, y ) ≤ 2T Y for all (t, y ) ∈ D.
∂y ∂y

so f (t, y ) = t y 2 satisfies Lipschitz condition with L = 2T Y .

What is going on with f (t, y ) = t y 2 ?
I f (t, y ) = t y 2 satisfies Lipschitz condition on the region
D = {(t, y ) | 0 ≤ t ≤ T, −Y ≤ y ≤ Y } .
I f (t, y ) = t y 2 doesn’t satisfy Lipschitz condition on region
D = {(t, y ) | 0 ≤ t ≤ T }.
Initial value problem
y 0 (t) = t y 2 (t), y (t0 ) = α > 0
has unique, but unbounded solution
2α
y (t) = ,
2 + α(t02 − t 2 )
the denominator of which vanishes at
r
2
t= + t02 .
α
What is going on with f (t, y ) = t y 2 ?
Initial value problem
y 0 (t) = t y 2 (t), y (t0 ) = α > 0
has unique, but unbounded solution
2α
y (t) = ,
2 + α(t02 − t 2 )
the denominator of which vanishes at
r
2
t= + t02 .
α
I for |t0 | < T , ODE has unique solution on
D = {(t, y ) | 0 ≤ t ≤ T, −Y ≤ y ≤ Y } .
q q
2
I for α + t02 < T ODE solution breaks down at t = α2 + t02
on
D = {(t, y ) | 0 ≤ t ≤ T }.