[PPT] - Lesson 2 Quadrature rules 1 Quadrature is the computation of PowerPoint Presentation

SLIDE 1

Lesson 2 Quadrature rules

1

SLIDE 2

The goal of standard quadrature methods is to express the integral as a finite sum

  so that the approximation converges:

Quadrature is the computation of integrals

 

Only in very special cases can integrals be computed exactly, otherwise, we must

compute them numerically (i.e., approximately)

Z b

a

f(x) dx

2

n

k=1

wn,kf(xn,k) lim

n→∞

n
i=1

wn,kf(xn,k) − b

a

f(x) dx

= 0

SLIDE 3

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

n = 3

The standard definition for an integral from

introductory calculus is  

SIMPLEST METHOD: RIGHTPOINT RULE

Z 1 cos x dx

f(h)

b

a

f(x) dx = lim

n→∞

b − a n

n

k=1

f

a + k b − a

n

h = b − a

n

SLIDE 4

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

n = 3

The standard definition for an integral from

introductory calculus is  

Idea: use the sum with finite n

This is equivalent to approximating the integral

by small rectangles

SIMPLEST METHOD: RIGHTPOINT RULE

Z 1 cos x dx

f(h)

b

a

f(x) dx = lim

n→∞

b − a n

n

k=1

f

a + k b − a

n

b

a

f(x) dx ≈ QR

n = h n

k=1

f(xk) for xk = a + k b − a n

h = b − a n

SLIDE 5

Alternatively, we can evaluate f on

the left:   

LEFTPOINT RULE

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

n = 3

Z 1 cos x dx

f(0)

b

a

f(x) dx ≈ QL

n = h n

k=1

f(xk) for xk = a + (k − 1)b − a n

h = b − a n

SLIDE 6

Instead of rectangles, we can use trapezoids
We then have

TRAPEZOIDAL RULE

Z 1 cos x dx

for and

6

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

n = 3

b

a

f(x) dx ≈

n−1

k=1

xk+1

xk

f(xk) + (x − xk)f(xk+1) − f(xk)

xk+1 − xk

dx
xk = a + (k − 1)h

h = b−a

n−1

4

SLIDE 7

Instead of rectangles, we can use trapezoids
We then have

TRAPEZOIDAL RULE

Z 1 cos x dx

for and

7

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

n = 3

b

a

f(x) dx ≈

n−1

k=1

xk+1

xk

f(xk) + (x − xk)f(xk+1) − f(xk)

xk+1 − xk

dx

=

n−1

k=1
f(xk)h + hf(xk+1) − f(xk)

2

n−1
xk = a + (k − 1)h

h = b−a

n−1

4

SLIDE 8

Instead of rectangles, we can use trapezoids
We then have

TRAPEZOIDAL RULE

Z 1 cos x dx

for and

8

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

n = 3

b

a

f(x) dx ≈

n−1

k=1

xk+1

xk

f(xk) + (x − xk)f(xk+1) − f(xk)

xk+1 − xk

dx

=

n−1

k=1
f(xk)h + hf(xk+1) − f(xk)

2

= h

2

n−1

k=1

[f(xk) + f(xk+1)]

xk = a + (k − 1)h h = b−a

n−1

4

SLIDE 9

TRAPEZOIDAL RULE

Z 1 cos x dx

9

0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

n = 3

Instead of rectangles, we can use trapezoids
We then have

for and xk = a + (k − 1)h h = b−a

n−1

b

a

f(x) dx ≈ QT

n = h

f(a)

2 +

n−1

k=2

f(xk) + f(b) 2

4

SLIDE 10

Experiments of quadrature error

10

SLIDE 11

We’ve introduced three quadrature rules

– Left point rule, right point rule and trapezoidal rule

Each use precisely n function evaluations
How fast will they converge to the true integral as n increases?
Let’s try some experiments

11

SLIDE 12

12

200 400 600 800 1000 n 0.1 0.2 0.3 0.4

Error approximating

1 x x = − 1

Right-point rule, Left-point rule, Trapezoidal rule

SLIDE 13

13

200 400 600 800 1000 n 0.1 0.2 0.3 0.4

Error approximating

1 x x = − 1

Right-point rule, Left-point rule, Trapezoidal rule

1000 2000 3000 4000 5000 n 10-6 10-4 0.01

Log-scale plot

SLIDE 14

14

200 400 600 800 1000 n 0.1 0.2 0.3 0.4

Error approximating

1 x x = − 1

Right-point rule, Left-point rule, Trapezoidal rule Log-log-scale plot

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

SLIDE 15

15

200 400 600 800 1000 n 0.1 0.2 0.3 0.4

Error approximating

1 x x = − 1

Right-point rule, Left-point rule, Trapezoidal rule Log-log-scale plot

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

1 n 1 n2

SLIDE 16

16

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x

SLIDE 17

16

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x 500x

100 200 500 1000 2000 5000 n 10-7 10-5 0.001 0.1

SLIDE 18

16

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x 500x

100 200 500 1000 2000 5000 n 10-7 10-5 0.001 0.1

|x − .1|

100 200 500 1000 2000 5000 n 10-7 10-5 0.001 0.1

SLIDE 19

16

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x 500x

100 200 500 1000 2000 5000 n 10-7 10-5 0.001 0.1 100 200 500 1000 2000 5000 n 10-7 10-5 0.001 0.1

(x − .1) |x − .1|

100 200 500 1000 2000 5000 n 10-7 10-5 0.001 0.1

SLIDE 20

17

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x

SLIDE 21

17

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x

5 10 20 50 100 n 10-15 10-12 10-9 10-6 0.001 1

10πx

SLIDE 22

17

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x

5 10 20 50 100 n 10-15 10-12 10-9 10-6 0.001 1

10πx 50πx

5 10 20 50 100 n 10-15 10-12 10-9 10-6 0.001 1

SLIDE 23

17

100 200 500 1000 2000 5000 n 10-6 10-4 0.01

x

5 10 20 50 100 n 10-15 10-12 10-9 10-6 0.001 1

10πx 2πx

5 10 20 50 100 n 10-13 10-10 10-7 10-4 0.1

50πx

5 10 20 50 100 n 10-15 10-12 10-9 10-6 0.001 1

SLIDE 24

Warning!

Only use these methods for periodic functions!!
Much better methods exist for non-periodic smooth functions
We will see this later in the course
The better method is not Simpson’s rule

18

SLIDE 25

Convergence of quadrature rules

19

SLIDE 26

20

We want to prove the observed behaviour

Left and right point rules converge like O(1/n) Trapezoidal rule converges like O

1/n2

All rules converge a lot faster for periodic and smooth functions

The key is integration by parts

SLIDE 27

Right-hand rule

21

Assume for simplicity that f is smooth (infinitely differentiable in [a, b])
We start with the error in one panel:

xk+1

xk

[f(x) − f(xk+1)] x = xk+1

xk

(x − xk) [f(x) − f(xk+1)] x

SLIDE 28

Right-hand rule

22

Assume for simplicity that f is smooth (infinitely differentiable in [a, b])
We start with the error in one panel:

xk+1

xk

[f(x) − f(xk+1)] x = xk+1

xk

(x − xk) [f(x) − f(xk+1)] x = [(x − xk)(f(x) − f(xk+1))]xk+1

xk

− xk+1

xk

(x − xk)f (x) x

SLIDE 29

Right-hand rule

23

Assume for simplicity that f is smooth (infinitely differentiable in [a, b])
We start with the error in one panel:

xk+1

xk

[f(x) − f(xk+1)] x = xk+1

xk

(x − xk) [f(x) − f(xk+1)] x = [(x − xk)(f(x) − f(xk+1))]xk+1

xk

− xk+1

xk

(x − xk)f (x) x = −1 2 xk+1

xk

(x − xk)2 f (x) x

SLIDE 30

Right-hand rule

24

Assume for simplicity that f is smooth (infinitely differentiable in [a, b])
We start with the error in one panel:

xk+1

xk

[f(x) − f(xk+1)] x = xk+1

xk

(x − xk) [f(x) − f(xk+1)] x = [(x − xk)(f(x) − f(xk+1))]xk+1

xk

− xk+1

xk

(x − xk)f (x) x = −1 2 xk+1

xk

(x − xk)2 f (x) x

= −1 2

(x − xk)2f (x)

xk+1

xk

+ 1 2 xk+1

xk

(x − xk)2f (x) x

2

x

SLIDE 31

Right-hand rule

25

Assume for simplicity that f is smooth (infinitely differentiable in [a, b])
We start with the error in one panel:

xk+1

xk

[f(x) − f(xk+1)] x = xk+1

xk

(x − xk) [f(x) − f(xk+1)] x = [(x − xk)(f(x) − f(xk+1))]xk+1

xk

− xk+1

xk

(x − xk)f (x) x = −1 2 xk+1

xk

(x − xk)2 f (x) x

= −1 2

(x − xk)2f (x)

xk+1

xk

+ 1 2 xk+1

xk

(x − xk)2f (x) x = −h2 2 f (xk+1) + 1 2 xk+1

xk

(x − xk)2f (x) x

SLIDE 32

26

We can bound the first term
h2

2 f (xk+1)

≤ h2

2

x |f (x)| =

b − a n 2

x |f (x)| = O

1 n2

we can also bound the second term
Summing over all panels we get

SLIDE 33

27

We can bound the first term
h2

2 f (xk+1)

≤ h2

2

x |f (x)| =

b − a n 2

x |f (x)| = O

1 n2

we can also bound the second term
xk+1

xk

(x − xk)2f (x) x

≤

x |f (x)| h3 = O

1 n3

Summing over all panels we get

SLIDE 34

28

We can bound the first term
h2

2 f (xk+1)

≤ h2

2

x |f (x)| =

b − a n 2

x |f (x)| = O

1 n2

we can also bound the second term
xk+1

xk

(x − xk)2f (x) x

≤

x |f (x)| h3 = O

1 n3

Summing over all panels we get

b

a

f(x) x − QR

n = n

k=1

xk+1

xk

[f(x) − f(xk+1)] x =

n

k=1

O 1 n2

= O

1 n

SLIDE 35

Trapezoidal rule

29

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x
where

is some point in

Summing over all the panels shows

error

SLIDE 36

Trapezoidal rule

30

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

= − xk+1

xk

(x − xk)

f (x) − f(xk+1) − f(xk)

h

x
where

is some point in

Summing over all the panels shows

error

SLIDE 37

Trapezoidal rule

31

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

= − xk+1

xk

(x − xk)

f (x) − f(xk+1) − f(xk)

h

x

= −h2

f (xk+1) − f(xk+1) − f(xk)

h

+ 1

2 xk+1

xk

(x − xk)2f (x) x x

where

is some point in

Summing over all the panels shows

error

SLIDE 38

Trapezoidal rule

32

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

= − xk+1

xk

(x − xk)

f (x) − f(xk+1) − f(xk)

h

x

= −h2

f (xk+1) − f(xk+1) − f(xk)

h

+ 1

2 xk+1

xk

(x − xk)2f (x) x = −h3f (χk) + 1 2 xk+1

xk

(x − xk)2f (x) x = O 1 n3

where χk is some point in [xk, xk+1]
Summing over all the panels shows O

1

n2

error

SLIDE 39

Euler–McLaurin Formula

33

SLIDE 40

34

We now want to show that the trapezoidal rule converges much faster for periodic

smooth functions

This involves the following modifications of the previous argument:

Replace x − xk by x − xk+xk+1

2

so the error terms cancel Integrate by parts again Replace (x − xk)2 by a suitable quadratic so the next term also cancels Integrate by parts again, and so on

SLIDE 41

35

Define the Bernoulli polynomials Bk(x) inductively by

B0(x) = 1 B

k(x) = kBk1(x) and

1

0 Bk(x) x = 1

These satisfy for
where

are called Bernoulli numbers

The first few are

SLIDE 42

36

Define the Bernoulli polynomials Bk(x) inductively by

B0(x) = 1 B

k(x) = kBk1(x) and

1

0 Bk(x) x = 1

These satisfy for k ≥ 2

Bk(1) = Bk(0) + 1 kBk1(x) x = Bk(0) = Bk where Bk are called Bernoulli numbers

The first few are

SLIDE 43

37

Define the Bernoulli polynomials Bk(x) inductively by

B0(x) = 1 B

k(x) = kBk1(x) and

1

0 Bk(x) x = 1

These satisfy for k ≥ 2

Bk(1) = Bk(0) + 1 kBk1(x) x = Bk(0) = Bk where Bk are called Bernoulli numbers

The first few are

B1(x) = x − 1 2

B1 := −1

2

B2(x) = x2 − x + 1

6

B2 = 1

6

B3(x) = x3 − 3

2x2 + x 2 (B3 = 0) B4(x) = x4 − 2x3 + x2 − 1 30

B4 = − 1

30

0.2

0.4 0.6 0.8 1.0

0.4
0.2

0.2 0.4 0.6 0.8 1.0

B0(x),B1(x),B2(x),B3(x),B4(x)

SLIDE 44

38

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

SLIDE 45

39

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

= −h xk+1

xk

B1 x − xk h f (x) − f(xk+1) − f(xk) h

x

SLIDE 46

40

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

= −h xk+1

xk

B1 x − xk h f (x) − f(xk+1) − f(xk) h

x

= −h2 2 B2 [f (xk+1) − f (xk)] + h2 2 xk+1

xk

B2 x − xk h

f (x) x

x

SLIDE 47

41

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

= −h xk+1

xk

B1 x − xk h f (x) − f(xk+1) − f(xk) h

x

= −h2 2 B2 [f (xk+1) − f (xk)] + h2 2 xk+1

xk

B2 x − xk h

f (x) x

= −h2 2 B2 [f (xk+1) − f (xk)] − h3 6 xk+1

xk

B3 x − xk h

f (x) x

SLIDE 48

42

Again start with the error in one panel:

xk+1

xk

f(x) − f(xk) − (x − xk)f(xk+1) − f(xk)

h

x

= −h xk+1

xk

B1 x − xk h f (x) − f(xk+1) − f(xk) h

x

= −h2 2 B2 [f (xk+1) − f (xk)] + h2 2 xk+1

xk

B2 x − xk h

f (x) x

= −h2 2 B2 [f (xk+1) − f (xk)] − h3 6 xk+1

xk

B3 x − xk h

f (x) x

. . . = −

m

k=2

hkBk k!

f (k1)(xk+1) − f (k1)(xk)
− hm

m! xk+1

xk

Bm x − xk h

f (m)(x) x

SLIDE 49

43

0.2 0.4 0.6 0.8 1.0 x

2.95
2.90
2.85
2.80
2.75

100 1000 500 200 2000 300 150 1500 700 n 10-11 10-8 10-5 0.01

ex(3x − x2 − 3)

1 n 1 n2 1 n4

SLIDE 50

43

0.2 0.4 0.6 0.8 1.0 x

2.95
2.90
2.85
2.80
2.75

100 1000 500 200 2000 300 150 1500 700 n 10-11 10-8 10-5 0.01

ex(3x − x2 − 3)

5 10 20 50 100 n 10-12 10-9 10-6 0.001 1

0.2 0.4 0.6 0.8 1.0 x 1.0 1.5 2.0 2.5 3.0 3.5

x +

x − 1

2 2 sinh

x − 1

2

+ ecos 2πx