Numerical Differentiation and Integration  Numerical Differentiation ◦ ◦ ◦ ◦

Finite Differences Interpolating Polynomials Taylor Series Expansion Richardson Extrapolation

 Numerical Integration ◦ Basic Numerical Integration ◦ Improved Numerical Integration ⇒ Trapezoidal, Simpson’s Rules ◦ Rhomberg Integration

ITCS 4133/5133: Numerical Comp. Methods

1

Numerical Differentiation and Integration

Numerical Differentiation and Integration  Many engineering applications require numerical estimates of derivatives of functions

 Especially true, when analytical solutions are not possible  Differentiation: Use finite differences  Integration (definite integrals): Weighted sum of function values at specified points (area under the curve).

ITCS 4133/5133: Numerical Comp. Methods

2

Numerical Differentiation and Integration

Application:Integral of a Normal Distribution 2

◦ Represented as a Gaussian, a scaled form of f (x) = e−x , very important function in statistics

◦ Not easy to determine indefinite integral - use numerical techniques Z

b

e−x

A=

2

a

ITCS 4133/5133: Numerical Comp. Methods

3

Numerical Differentiation and Integration

Application:Integral of a Sinc f (x) =

ITCS 4133/5133: Numerical Comp. Methods

sin(x) x

4

Numerical Differentiation and Integration

Numerical Differentiation:Approach ⇒ Evaluate the function at two consecutive points, separated by ∆x, of the independent variable, and taking their difference:

df (x) f (x + ∆x) − f (x) ≈ dx ∆x ⇒ Fit a function to a set of points that define the relationship between the independent and dependent variables, such as an nth order polynomial; then differentiate the polynomial

ITCS 4133/5133: Numerical Comp. Methods

5

Numerical Differentiation and Integration

Finite Differences Forward Difference df (x) f (x + ∆x) − f (x) = dx ∆x Backward Difference df (x) f (x) − f (x − ∆x) = dx ∆x Two Step Method, Central Difference df (x) f (x + ∆x) − f (x − ∆x) = dx 2∆x ITCS 4133/5133: Numerical Comp. Methods

6

Numerical Differentiation and Integration

Finite Differences: Example

ITCS 4133/5133: Numerical Comp. Methods

7

Numerical Differentiation and Integration

Forward Differences: Derivation  Finite difference formulas can be derived from the Taylor series. h2 00 f (x + h) = f (x) + hf (x) + f (η), 2 0

For h = xi+1 − xi ,

f 0(xi) =

f (xi+1) − f (xi) h 00 − f (η), h 2

where xi < η < xi+1 .

ITCS 4133/5133: Numerical Comp. Methods

8

Numerical Differentiation and Integration

Backward Differences: Derivation h2 00 f (x + h) = f (x) + hf (x) + f f (η), 2 0

For h = xi−1 − xi ,

h2 00 f (xi + h) = f (xi + xi−1 − xi) = f (xi−1) = f (xi) + hf (xi) + f f (η) 2 f (xi−1) − f (xi) h 00 f 0(xi) = − f (η) h 2 0

where xi−1 < η < xi .

ITCS 4133/5133: Numerical Comp. Methods

9

Numerical Differentiation and Integration

Central Differences: Derivation Use the next higher order Taylor polynomial,

h2 00 h3 000 f (xi+1) = f (xi) + hf (xi) + f (xi) + f (η1), 2 6 2 h 00 h3 000 0 f (xi−1) = f (xi) − hf (xi) + f (xi) − f (η2) 2 6 0

with x < η1 < x + h, x − h < η2 < x Thus,

h3 000 f (xi+1) − f (xi−1) = 2hf (xi) + [f (η1) + f 000(η2)] 6 0

or,

f (xi+1) − f (xi−1) h2 000 f (xi) = − f (η), xi−1 < η < xi + 1 2h 6 0

ITCS 4133/5133: Numerical Comp. Methods

10

Numerical Differentiation and Integration

Second Derivatives h3 000 h4 (4) h2 00 f (x + h) = f (x) + hf (x) + f (x) + f (x) + f (x)(η1) 2 3! 4! 2 3 h 00 h 000 h4 (4) 0 f (x − h) = f (x) − hf (x) + f (x) − f (x) + f (x)(η2) 2 3! 4! 0

Thus

h4 (4) f (x + h) + f (x − h) = 2f (x) + h f (x) + [f (x)(η1) + f (4)(x)(η2)] 4! 1 f 00(x) ≈ 2 [f (x + h) − 2f (x) + f (x − h)] h 2 00

with truncation error of the O(h4 ) ITCS 4133/5133: Numerical Comp. Methods

11

Numerical Differentiation and Integration

Partial Derivatives  Generally, interested in partial derivatives of functions of 2 variables, (xi, yj ), based on a mesh of points.  Subscripts denote partial derivatives.

ITCS 4133/5133: Numerical Comp. Methods

12

Numerical Differentiation and Integration

Partial Derivatives (contd)  Laplacian operator: ∆2u = uxx + uyy  Biharmonic Operator: ∆4u = uxxxx + uxxyy + uyyyy

ITCS 4133/5133: Numerical Comp. Methods

13

Numerical Differentiation and Integration

Using Interpolating Polynomials  Given a function in discrete form, the data can be interpolated to fit an nth order polynomial

 For those functions that are in analytical form, but difficult to differentiate, the function can be discretized and fitted by a polynomial.

f (x) = bnxn + bn−1xn−1 + · · · + b1x + b0 whose derivative is

df (x) = nbnxn−1 + (n − 1)bn−2xn−2 + · · · + b1 dx

ITCS 4133/5133: Numerical Comp. Methods

14

Numerical Differentiation and Integration

Power Series Type Interpolating Polynomial Consider the nth order polynomial passing through (n+1) points; the (n+ 1) coefficients can be uniquely determined.

Example f (x) = a0 + a1x + a2x2 Consider f (xi ), f (xi + h), f (xi + 2h)

f (xi) = a0 + a1xi + a2x2i 2 f (xi + h) = a0 + a1(xi + h) + a2(xi + h) 2 f (xi + 2h) = a0 + a1(xi + 2h) + a2(xi + 2h)

ITCS 4133/5133: Numerical Comp. Methods

15

Numerical Differentiation and Integration

Example (contd) For the 3 points, xi = 0, h, 2h,

fi = a0 fi+1 = a0 + a1h + a2h2 fi+2 = a0 + 2a1h + 4a2h2 which has the solution,

a0 = f i −fi+2 + 4fi+1 − 3fi a1 = 2h fi+2 − 2fi+1 + fi a2 = 2h2

ITCS 4133/5133: Numerical Comp. Methods

16

Numerical Differentiation and Integration

Example (contd) f (x) = a0 + a1x + a2x2 0

f 0(xi) = fi = f 0(xi = 0) = a1 + 2a2xi = a1 −fi+2 + 4fi+1 − 3fi = 2h Similarly, second derivatives can also be obtained.

ITCS 4133/5133: Numerical Comp. Methods

17

Numerical Differentiation and Integration

Numerical Integration  Integration can be thought of as considering some continuous function f (x) and the area A subtended by it; for instance, within a particular interval b

Z

f (x)dx

A= a

 Numerical Integration is needed when f (x) does not have a known analytical solution, or, if f(x) is only defined at discrete points.

 There are two approaches to a numerical solution: ⇒ Fitting polynomials to f (x) and integrating using analytical calculus, ⇒ Area under f (x) can be approximated using geometric shapes defined by adjacent points. ITCS 4133/5133: Numerical Comp. Methods

18

Numerical Differentiation and Integration

Interpolation Approach We can derive an nth order polynomial (for instance, Gregory-Newton method), which has the general form as

f (x) = b1xn + b2xn−1 + · · · + bnx + bn+1 After the bi s are determined, this can be integrated as

Z

b1xn+1 b2xn bn x 2 f (x)dx = + + ··· + + bn+1x n+1 n 2

which can be solved for, within the limits of the integration

ITCS 4133/5133: Numerical Comp. Methods

19

Numerical Differentiation and Integration

Basic Numerical Integration (Newton-Cotes Closed Formulas)  Trapezoid Rule: Approximates function by a straightline to compute area under curve.

Z a

b

h f (x)dx ≈ [f (x0) + f (x1)] 2

 Solution is exact, for polynomials of degree ≤ 1, i.e. linear functions.

ITCS 4133/5133: Numerical Comp. Methods

20

Numerical Differentiation and Integration

Basic Numerical Integration (Newton-Cotes Closed Formulas)  Simpson’s Rule: Instead of a linear approximation, use a quadratic polynomial:

h=

b−a b+a , x0 = a, x1 = x0 + h = , x2 = b 2 2

Approximate integral is given by

Z a

b

h f (x)dx ≈ [f (x0) + 4f (x1) + f (x2) 3

ITCS 4133/5133: Numerical Comp. Methods

21

Numerical Differentiation and Integration

Basic Simpson’s Rule: Examp1e 1

ITCS 4133/5133: Numerical Comp. Methods

22

Numerical Differentiation and Integration

Basic Simpson’s Rule: Examp1e 2

ITCS 4133/5133: Numerical Comp. Methods

23

Numerical Differentiation and Integration

Improved Numerical Integration  Improve the accuracy by applying lower order methods repeatedly on several subintervals.

 Known as composite integration  Simpson, Trapezoid rules use equal subintervals; using unequal (adaptive) intervals leads to Gaussian Quadrature methods.

ITCS 4133/5133: Numerical Comp. Methods

24

Numerical Differentiation and Integration

Composite Trapezoid Rule  Approximates the area under the function f (x) by a set of discrete trapezoids fitted between each pair of points of the dependent variable (and the X axis)

Z

xn

x1

n−1 X f (xi+1 + f (xi) f (x)dx ≈ (xi+1 − xi) 2 i=1

f(x)

f(x)

x x1

x2

x3

x4

If the intervals are the same, h = (b − a)/n, we get

Z

xn

f (x)dx ≈ x1

b−a [f (a) + 2f (x1) + . . . + f (b)] 2n

ITCS 4133/5133: Numerical Comp. Methods

25

Numerical Differentiation and Integration

Composite Trapezoid Rule:Example

ITCS 4133/5133: Numerical Comp. Methods

26

Numerical Differentiation and Integration

Composite Trapezoid Rule:Algorithm

ITCS 4133/5133: Numerical Comp. Methods

27

Numerical Differentiation and Integration

Composite Trapezoid Rule: Notes ⇒

f (xi+1 +f (xi ) 2

represents the average height of the trapezoid.

⇒ Linear approximation between pairs of successive points used; error is proportional to the distance between sample points.

Error ⇒ Error in the trapezoidal can be approximated by Taylor’s series xi+1 − xi 00 f (x) 2 where f 00 (x) is evaluated at the value of x that maximizes the second derivative between the two sample points.

ITCS 4133/5133: Numerical Comp. Methods

28

Numerical Differentiation and Integration

Composite Simpson’s Rule  Uses the same idea as the composite Trapezoid rule, by using multiple sub-intervals to improve the Simpson rule approximation.

 For 2 subintervals, Simpson’s rule has the following form: Z b Z x2 Z b f (x)dx = f (x)dx + f (x)dx a

a

x2

h h [f (a) + 4f (x1) + f (x2)] + [f (x2) + 4f (x3) + f (b)] 3 3 h ≈ [f (a) + 4f (x1) + 2f (x2) + 4f (x3) + f (b)]) 3

≈

Trapezoidal Rule

f(x)

f(x) Simpson’s Rule

ITCS 4133/5133: Numerical Comp. Methods a

29

(a+b)/2

x Numerical Differentiation and Integration b

Composite Simpson’s Rule (contd)  For n (n must be even) subintervals, h = (b − a)/n we get Z b h f (x)dx = [f (a) + 4f (x1) + 2f (x2) + 4f (x3) + 3 a 2f (x4) + . . . + 2f (xn−2 + 4f (xn−1) + f (b)  n must be even number of subintervals.  Error Bound: Based on the 4th derivative: 5

(xi+1 − xi) (4) f (xi) Error = 90 where f (4) (xi ) is evaluated at the value of x that maximizes the 4th derivative between the two sample points. ITCS 4133/5133: Numerical Comp. Methods

30

Numerical Differentiation and Integration

Composite Simpson’s Rule:Example

ITCS 4133/5133: Numerical Comp. Methods

31

Numerical Differentiation and Integration

Composite Simpson’s Rule:Algorithm

ITCS 4133/5133: Numerical Comp. Methods

32

Numerical Differentiation and Integration

Romberg Integration  Simpson’s rule improves on Trapezoid rule with additional evaluations of f(x), required to fit a more accurate polynomial

 More evaluations will further improve the integral, leading to the Romberg Integration

I01 =

b−a (f (a) + f (b)) 2

A second estimate can be obtained with subdividing the interval ab into 2 equal intervals is given by

I11

  1 b−a 1 = f (a) + f (m) + f (b) 2 2 2    1 b−a = I01 + (b − a)f a + 2 2

ITCS 4133/5133: Numerical Comp. Methods

33

Numerical Differentiation and Integration

Romberg Integration(contd.) A third estimate, I21 , using 3 intermediate points m1 , m2 , m3 is given by

I21

  b−a 1 1 1 1 1 = f (a) + f (m1) + f (m2) + f (m3) + f (b) 2 4 2 2 2 4 " #   3 1 b−a X b−a = I11 + f a+ k 2 2 k=1,k6=2 4

This leads to the recursive relationship

 # 2i −1 X 1 b−a b−a Ii1 = Ii−1,1 + i−1 f a+ i k , for i = 1, 2, . . . 2 2 2 k=1,3,5 "

ITCS 4133/5133: Numerical Comp. Methods

34

Numerical Differentiation and Integration

Romberg Integration(contd.)  # 2i −1 X 1 b−a b−a Ii1 = Ii−1,1 + i−1 f a+ i k , for i = 1, 2, . . . 2 2 2 k=1,3,5 "

We can combine these estimates using Richardson’s formula,

4j−1Ii+1,j−1 − Ii,j−1 Iij = 4j−1 − 1 for i = 0, 1, . . . , N − j + 1, j = 1, 2, . . . N . The estimates form an upper triangular matrix:

I

01

 II11  21  I31  .   . .

I02 I12 I22 . . .

I03 I13 . . .

I04 . . . .

··· ··· ··· ··· ···

IN −2,3

IN −1,2



I0,N −1 I0,N I0,N +1 I1,N −1 I1,N  I2,N −1 

   

IN 1 ITCS 4133/5133: Numerical Comp. Methods

35

Numerical Differentiation and Integration

Romberg Integration(contd.) Algorithm (To solve for I0N  Determine I01 (Trapezoid formula)  Determine Ii1, for i = 1, 2, . . .  Determine Ii2, Ii3, etc., using Richardson’s formula.

ITCS 4133/5133: Numerical Comp. Methods

36

Numerical Differentiation and Integration

ITCS 4133/5133: Numerical Comp. Methods

37

Numerical Differentiation and Integration