MATLAB

1

MIE 124

Numerical Integration

I – Overview

1 – Intro

 Very often in engineering we have a function or data and want to know the integral of that function.

 Examples include:

o Measure forces along a wind turbine blade, and want to compute the bending moment on the root.

o Measure the flow velocity at many points in a river, and want to calculate the total volume or mass flow

rate.

o Measure the heat addition and heat loss to a building over time, and want to calculate the total heat flux.

2 – Math refresher

 We have a function f(x), and we want to determine the integral of that function between limits a and b.

( ) ( ) ( ) b

b

I f x dx F b F a  

 F(x) is the integral of f(x), and I is the value of the integral between limits a and b.

 \\ Integration is related to summation, that’s why the integral sign is a fancy S.

 Graphically, I is the area underneath the curve f(x).

 [] Draw F

 You probably all remember how to integrate a polynomial:

2

1

1

bb n n

a a

x x dx

n



 

 But sometimes we have function where we can’t determine the integral.

 Numerical integration is the solution.

2 – Numerical Integration intro

 Say we have a complicated function that cannot be integrated analytically.

 [] Draw (a)

 In Num. Int, we convert the continuous function to discrete values of the function at a number of points.

 [] Draw (b)

 Then we combine the discrete values using some sort of weighting system to estimate the integral.

 [] Draw (c)

3 – Graphical method

 An alternative is a graphical method.

 [] Draw F

3

 Tedious and inexact.

II – Trapezoid Method

1 – Graphically

 \\ Probably many of you are familiar with the trapezoid method.

 Graphically, we connect points on f(x) using straight-line segments, then sum the area in the resulting

trapezoids.

 [] Draw F

2 – Mathematically

 For a single segment:

4

( ) ( ) ( )

2

f a f b I b a

  

 For multiple segments, let’s assume n+1 equally spaced base points (x1, x2, …, xn+1), and n equal segments.

 The segment width is h=(b-a)/n.

 \\ Take 2 minutes, apply the trapezoid rule for each segment, and group terms. Try to determine an expression

that uses a summation over an arbitrary number of segments.

 Results is:

1 1 2

( ) 2 ( ) ( ) 2

n

i n i

h I f x f x f x 



     

  

 Note that the reading I posted uses indices from 0-n, but since first index in Matlab is 1, this is easier.

3 – Error analysis

 Can derive expression for the error when using the trapezoid rule.

 Result is E α 1/n^2

 If double n, the error decreases by a factor of 4.

 If the function is linear, what is the error?

 Exact for a linear function, trapezoid rule is a “first order” solution.

 How do you know if n is big enough?

 Can check approximate error, i.e. absolute relative change in I as n increases (just like root finding).

4 – Example and Numerical solution

 For an example, let’s use the expression for the velocity of a falling object as a function of time:

 ( / )( ) 1 c m tgmv t e c

 

 We want to calculate the distance D that the object falls between t1=0 s and tf=50 s.

 [] First go to Matlab and plot

 \\ Now take 5 minutes and write pseudocode to integrate v(t) numerically to find D, using n=16.

 [] Go to Matlab and show code, then run code and calculate D.

 I used a subfunction to calculate v(t), more convenient than entering the equation every time.

 Linspace is useful to generate equally spaced segments between a and b.

5 – Example and effect of discretization

 Let’s investigate the effect of increasing n.

 Use n = [2^4; 2^8; 2^12; 2^16; 2^20] (see effect of increasing by factors of 2).

5

 [] Open num_int_ex_moodle in Matlab and walk through code

 Note, using a loop to store values of integral for each value of n.

 [] Discuss plots.

 \\ You will do something similar for your homework except using Simpson’s method

III – Simpson’s Method

 Instead of using more segments in trapezoid rule to get a better estimate of I, use a higher order method.

 Simpson’s method connects 3 points using a parabola.

1 – Graphically

 Graphically, we connect points on f(x) using parabolas, then sum the area.

 [] Draw F

2 – Mathematically

 For a single segment:

( ) 4 f((a b) / 2) ( ) ( )

6

f a f b I b a

    

 Again assume n+1 equally spaced base points (x1, x2, …, xn+1), n equal segments, and h=(b-a)/n.

 Multiple segments:

1

1 1 2,4,6 3,5,7

( ) 4 ( ) 2 ( ) ( ) 3

n n

i i n i i

h I f x f x f x f x



  

      

   

 Assume n is even (could use trapezoid method on first segment if odd).

 Again reading uses 0-n, but since first index in Matlab is 1, this is easier.

3 – Error analysis

 Can derive expression for the error when using the Simpsons.

 Result is E α 1/n^4

6

 If double n, the error decreases by a factor of 16.

 If the function is quadratic, is there error? No.

 Actually it is exact for a linear, quadratic, or cubic functions, Simpsons rule is a “second order” solution but

“third order” accurate.

4 – Example and Numerical solution

 Use same example:

 ( / )( ) 1 c m tgmv t e c

 

 Compare Trapezoid and Simpsons.

 \\ Won’t show you Simpsons code, almost identical to trapezoid code. You will implement for your

homework.

 [] Run num_int_ex.

 [] Discuss plots.

5 – In class exercise on separate document