We have seen two applications: signal smoothing root finding - - PowerPoint PPT Presentation

DIFFERENTIATION AND INTEGRATION

Lesson 10

• We have seen two applications:

– signal smoothing – root finding

• Today we look

– differentation – integration

• These will form the basis for solving ODEs

DIFFERENTIATION OF FOURIER SERIES

• For functions on the periodic interval, we have the Fourier representation

f(θ) =

• k=

ˆ fkkθ

• Its derivative is (formally) obvious:

f (θ) ∼

• k=

k ˆ fkkθ

• When does this converge? Whenever ˆ

fk decays fast enough

• Numerically, we obtain the approximation:
• ...
• where

and is the DFT

• For functions on the periodic interval, we have the Fourier representation

f(θ) =

• k=

ˆ fkkθ

• Its derivative is (formally) obvious:

f (θ) ∼

• k=

k ˆ fkkθ

• When does this converge? Whenever ˆ

fk decays fast enough

• Numerically, we obtain the approximation:

f (θ) ≈

• αθ, . . . , βθ
• α

... β

• Ff

where f = (f(θ1), . . . , f(θm)) and F is the DFT

Pointwise convergence of derivative of at zero

500 1000 1500 2000 10-12 10-9 10-6 0.001 1

Numerical derivative Direct interpolation number of points number of points

f(θ) = ecos(10θ−1)

500 1000 1500 2000 10-13 10-10 10-7 10-4 0.1

• The Nth order derivative is:

f (N)(θ) ∼

• k=

(k)N ˆ fkkθ

• Numerically, we obtain the approximation:

f (N)(θ) ≈

• αθ, . . . , βθ
• α

... β

• N

Ff where f = (f(θ1), . . . , f(θm)) and F is the DFT

Pointwise convergence of 10th derivative of at zero Numerical 10th derivative Direct interpolation number of points number of points

f(θ) = ecos(10θ−1)

500 1000 1500 2000 10-8 10-6 10-4 0.01 1 200 400 600 800 1000 10-14 10-11 10-8 10-5 0.01 10

INTEGRATION OF FOURIER SERIES

• A function define don the periodic interval has the indefinite integral
• f θ =
• k=,k=0

ˆ fk k kθ + ˆ f0θ + C

• This will converge whenever the Fourier series does!
• Numerically, we obtain the approximation:
• ...
• ...
• where

and is the DFT

• This is stable because the error in each computed

is multiplied by a bounded number

• A function define don the periodic interval has the indefinite integral
• f θ =
• k=,k=0

ˆ fk k kθ + ˆ f0θ + C

• This will converge whenever the Fourier series does!
• Numerically, we obtain the approximation:
• f θ ≈
• αθ, . . . , βθ
• 1

α

...

1

• 1
• ...

1 β

• Ff +θe

0 Ff +C

where f = (f(θ1), . . . , f(θm)) and F is the DFT

• This is stable because the error in each computed ˆ

fk is multiplied by a bounded number

Pointwise convergence of integral of at zero number of points

200 400 600 800 1000 10-15 10-11 10-7 0.001

f(θ) = ecos(10θ−1)

• Letting C = 0 at each stage, we can iterate this N times:
• · · ·
• f θN ≈
• αθ, . . . , βθ
• 1

α

...

1 − 1

• ...

1 β

• N

Ff+ ˆ f0 N!θN

• The stability of this approximation is maintained

DIFFERENTIATION OF TAYLOR SERIES

• For functions in the disk, we have the Taylor series

f(z) =

• k=0

ˆ fkzk

• Derivative :

f (z) =

• k=0

k ˆ fkzk1

• Numerically, we obtain the (m − 1) × m matrix approximation:

f (z) ≈

• 1 | · · · | zm2
• 1

2 ... m − 1

• T f

where f = (f(z1), . . . , f(zm)) and T is the discrete Taylor transform

50 100 150 200 10-14 10-11 10-8 10-5 0.01

First derivative

50 100 150 200 10-8 10-5 0.01 10 104 107

10th derivative Error approximating exp(z) for z = {.1,.5,1.}exp(.1i) number of points

INTEGRATION OF TAYLOR SERIES

• Integral :
• f(z) z =
• k=0

ˆ fk k + 1zk+1 + C

• Numerically, we obtain the (m + 1) × m matrix approximation:
• f(z) z ≈

1 | · · · | zm

• 1

1 2

...

1 m+1

• T f + C

where f = (f(z1), . . . , f(zm)) and T is the discrete Taylor transform

• Now this is stable both on and inside the unit circle!

First integral 10th integral Error approximating exp(z) for z = {.1,.5,1.}exp(.1i) number of points

50 100 150 200 10-14 10-11 10-8 10-5 0.01 50 100 150 200 10-16 10-14 10-12 10-10 10-8

DIFFERENTIATION OF LAURENT SERIES

• For functions on the circle, we have the Laurent series

f(z) =

• k=

ˆ fkzk

• Derivative :

f (z) =

• k=

k ˆ fkzk1 =

• k=,k=1

(k + 1) ˆ fk+1zk

• Numerically, we obtain the (square) approximation, where we changed our basis:

f (z) ≈

• zα1 | · · · | zβ1
• α

... β

• Ff

where f = (f(z1), . . . , f(zm))

• Clearly, it will only be somewhat accurate on the unit circle, and the error will grow

INTEGRATION OF LAURENT SERIES

• For functions on the circle, we have the Laurent series

f(z) =

• k=

ˆ fkzk

• Integral:
• f(z) z =
• k=,k=1

1 k + 1 ˆ fkzk+1 + ˆ f1 z + C

• Numerically, we obtain the approximation where we changed bases:
• f(z) z ≈
• zα+1 | · · · | zβ+1
• 1

α+1

...

1 1

1 ...

1 β+1

• Ff+e1Ff z+C

where f = (f(z1), . . . , f(zm))

• This will be stable
• However, we cannot (easily) iterate integrals!

INTEGRATION OF CHEBYSHEV SERIES

• We want to compute
• f(x) x =
• k=0

ˇ fk

• Tk(x) x
• The first two moments are
• T0(x) x = x+C = T1(x)+C

and

• T1(x) x = x2

2 +C = T2(x) − T0(x) 4 +C

• For

, we do the change of variables to map to the unit circle:

• (Now

!)

• We want to compute
• f(x) x =
• k=0

ˇ fk

• Tk(x) x
• The first two moments are
• T0(x) x = x+C = T1(x)+C

and

• T1(x) x = x2

2 +C = T2(x) − T0(x) 4 +C

• For k > 1, we do the change of variables x = J(z) to map to the unit circle:

x

a

Tk(x) x = J−1

↓

(x) J−1

↓

(a)

Tk(J(z))J(z) z = 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk

1 − 1 z2

• z
• (Now

!)

• We want to compute
• f(x) x =
• k=0

ˇ fk

• Tk(x) x
• The first two moments are
• T0(x) x = x+C = T1(x)+C

and

• T1(x) x = x2

2 +C = T2(x) − T0(x) 4 +C

• For k > 1, we do the change of variables x = J(z) to map to the unit circle:

x

a

Tk(x) x = J−1

↓

(x) J−1

↓

(a)

Tk(J(z))J(z) z = 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk

1 − 1 z2

• z

= 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk − zk2 − zk2

z

• (Now

!)

• We want to compute
• f(x) x =
• k=0

ˇ fk

• Tk(x) x
• The first two moments are
• T0(x) x = x+C = T1(x)+C

and

• T1(x) x = x2

2 +C = T2(x) − T0(x) 4 +C

• For k > 1, we do the change of variables x = J(z) to map to the unit circle:

x

a

Tk(x) x = J−1

↓

(x) J−1

↓

(a)

Tk(J(z))J(z) z = 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk

1 − 1 z2

• z

= 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk − zk2 − zk2

z = 1 4 zk+1 k + 1 + z1k 1 − k − zk1 k − 1 − zk1 −k − 1

• + C

(Now z = J1

(x)!)

• We want to compute
• f(x) x =
• k=0

ˇ fk

• Tk(x) x
• The first two moments are
• T0(x) x = x+C = T1(x)+C

and

• T1(x) x = x2

2 +C = T2(x) − T0(x) 4 +C

• For k > 1, we do the change of variables x = J(z) to map to the unit circle:

x

a

Tk(x) x = J−1

↓

(x) J−1

↓

(a)

Tk(J(z))J(z) z = 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk

1 − 1 z2

• z

= 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk − zk2 − zk2

z = 1 4 zk+1 k + 1 + z1k 1 − k − zk1 k − 1 − zk1 −k − 1

• + C

(Now z = J1

(x)!)

= 1 2(k + 1) zk+1 + zk1 2 − 1 2(k − 1) zk1 + z1k 2 + C

• We want to compute
• f(x) x =
• k=0

ˇ fk

• Tk(x) x
• The first two moments are
• T0(x) x = x+C = T1(x)+C

and

• T1(x) x = x2

2 +C = T2(x) − T0(x) 4 +C

• For k > 1, we do the change of variables x = J(z) to map to the unit circle:

x

a

Tk(x) x = J−1

↓

(x) J−1

↓

(a)

Tk(J(z))J(z) z = 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk

1 − 1 z2

• z

= 1 4 J−1

↓

(x) J−1

↓

(a)

• zk + zk − zk2 − zk2

z = 1 4 zk+1 k + 1 + z1k 1 − k − zk1 k − 1 − zk1 −k − 1

• + C

(Now z = J1

(x)!)

= 1 2(k + 1) zk+1 + zk1 2 − 1 2(k − 1) zk1 + z1k 2 + C = Tk+1(x) 2(k + 1) − Tk1(x) 2(k − 1) + C

• We thus have the integration formula
• f(x) x =
• k=0

ˇ fk

• Tk(x) x

= C + ˇ f0T1(x) + ˇ f1 T2(x) − T0(x) 4 + 1 2

• k=2

ˇ fk Tk+1(x) k + 1 − Tk1(x) k − 1

• Numerically, we approximate

and have the matrix for integration:

• ...

... where and is the discrete cosine transform (DCT)

• We thus have the integration formula
• f(x) x =
• k=0

ˇ fk

• Tk(x) x

= C + ˇ f0T1(x) + ˇ f1 T2(x) − T0(x) 4 + 1 2

• k=2

ˇ fk Tk+1(x) k + 1 − Tk1(x) k − 1

• = C −

ˇ f1 4 +

• ˇ

f0 − ˇ f2 2

• T1(x) + 1

2

• k=2

ˇ fk1 − ˇ fk+1 k

• Tk(x)
• Numerically, we approximate

and have the matrix for integration:

• ...

... where and is the discrete cosine transform (DCT)

• We thus have the integration formula
• f(x) x =
• k=0

ˇ fk

• Tk(x) x

= C + ˇ f0T1(x) + ˇ f1 T2(x) − T0(x) 4 + 1 2

• k=2

ˇ fk Tk+1(x) k + 1 − Tk1(x) k − 1

• = C −

ˇ f1 4 +

• ˇ

f0 − ˇ f2 2

• T1(x) + 1

2

• k=2

ˇ fk1 − ˇ fk+1 k

• Tk(x)
• Numerically, we approximate

f(x) ≈

• 1 | · · · | Tn1(x)
• Cf

and have the (n + 1) × n matrix for integration:

• f(x) x ≈
• 1 | · · · | Tn(x)
• − 1

4

1 − 1

2 1 4

− 1

4

... ...

1 2(n2)

−

1 2(n2) 1 2(n1) 1 2n

• Cf

where f = (f(x1), . . . , f(xn)) and C is the discrete cosine transform (DCT)

First integral 10th integral Error approximating exp(x) for x = .1 number of points

10 20 30 40 50 10-13 10-10 10-7 10-4 0.1 10 20 30 40 50 10-14 10-12 10-10 10-8

DIFFERENTIATION OF CHEBYSHEV SERIES

• Let's try decomposing the terms:

f (x) =

• k=0

ˇ fkT

k(x)

• We want to rewrite T

k(x) in terms of T0(x), . . . , Tk1(x)

• Let's try an experiment:

Use Tk(x) = k x, so T

k(x) = k k x

• 1x2

Therefore, we can numerically evaluate CTk(x) = C

• T

k(x1)

. . . T

k(xn)

• to get the coefficients

0 1. 3. 5. 7. 9. 4. 8. 12. 16. 20. 6. 10. 14. 18. 8. 12. 16. 20. 10. 14. 18. 12. 16. 20. 14. 18. 16. 20. 18. 20.

Problem: the operation is dense!

CT

0(x) | · · · | CT 10(x)

=

• Instead, we will use the fact that differentiation is the opposite of integration
• We want to find the vector of coefficients u = (ˇ

u0, . . . , ˇ un2) so that             − 1

4

1 − 1

2 1 4

− 1

4

... ...

1 2(n3)

−

1 2(n3) 1 2(n2) 1 2(n1)

            u = Cf

• We can apply backward substitution:

. . .

• What about the last condition

? This condition is not necessary because the constant of integration is arbitrary

• Instead, we will use the fact that differentiation is the opposite of integration
• We want to find the vector of coefficients u = (ˇ

u0, . . . , ˇ un2) so that             − 1

4

1 − 1

2 1 4

− 1

4

... ...

1 2(n3)

−

1 2(n3) 1 2(n2) 1 2(n1)

            u = Cf

• We can apply backward substitution:

ˇ un2 = 2(n − 1) ˇ fn1 ˇ un3 = 2(n − 2) ˇ fn2 . . .

• What about the last condition

? This condition is not necessary because the constant of integration is arbitrary

• Instead, we will use the fact that differentiation is the opposite of integration
• We want to find the vector of coefficients u = (ˇ

u0, . . . , ˇ un2) so that             − 1

4

1 − 1

2 1 4

− 1

4

... ...

1 2(n3)

−

1 2(n3) 1 2(n2) 1 2(n1)

            u = Cf

• We can apply backward substitution:

ˇ un2 = 2(n − 1) ˇ fn1 ˇ un3 = 2(n − 2) ˇ fn2 ˇ un4 = 2(n − 3) ˇ fn3 + ˇ un2 . . . ˇ u0 = ˇ f1 + ˇ u2 2

• What about the last condition − ˇ

u1 4 = ˇ

f0? This condition is not necessary because the constant of integration is arbitrary

First derivative 10th derivative Error approximating exp(x) for x = .1 number of points

10 20 30 40 50 10-13 10-10 10-7 10-4 0.1 10 20 30 40 50 10-5 10-4 0.001 0.01 0.1 1 10