numerical modeling.., Lecture notes of Numerical Methods in Engineering

Download numerical modeling.. and more Lecture notes Numerical Methods in Engineering in PDF only on Docsity!

Numerical Modeling in

Biomedical Systems

BME 125:

Lecture 10

Newton- Cotes

Differentiation

Interpolation

Finite difference methods

Finite Difference Methods, Interpolation, and Integration

Spline interpolation

Newton’s polynomials

Errors in FD methods

Integration

Simpson’s rules

-5 1 2 3 4 5 6 7 8

0

5

10

15

20

x

y

7th order polynomial

x = [1 2 3 4 5 6 7 8]; y = [4 7 5 3 8 0 3 2 ];

Example: 7th^ order polynomial fit to 8 data points

f  x a 1 x^2 a 2 xa 3 Careful – descending powers!

polyfit returns the coefficients p 1 , p 2 ,… in vector “p”. Function polyval takes

these coefficients along with a set of x values (xi) and calculates the polynomial f (x)

xi = [1:0.1:8]; yii = polyval(p,xi);

figure; plot(x,y,'ro',xi,yii); xlabel('x'); ylabel('y'); title('7th order polynomial');

p = polyfit(x,y,7);

20000 2002 2004 2006 2008 2010 2012

20

40

60

80

100

x

f(x)

Shifted fit

20000 2002 2004 2006 2008 2010 2012

20

40

60

80

100

x

f(x)

Raw fit - Large x-values

figure; plot(x,y,'ro',xi,yi); xlim([2000 2012]); ylim([0 100]); xlabel('x'); ylabel('f(x)'); title('Raw fit - Large x-values');

x = [2001 2002 2003 2004 2005 2006 2007 2009 2010 2011]; xi = x-2006; y = [18 28 23 28 26 24 50 28 36 39];

Shifting Data for Polynomial Fitting

(See note on Canvas > Resources > Supplementary Materials for an explanation of why polyfit can experience problems fitting high order polynomials to large x-values)

x = [2001 2002 2003 2004 2005 2006 2007 2009 2010 2011]; y = [18 28 23 28 26 24 50 28 36 39]; p = polyfit(x,y,9); xi = [2001:0.1:2011]; yi = polyval(p,xi); figure; plot(x,y,'ro',xi+2006,yi); xlim([2000 2012]); ylim([0 100]); xlabel('x'); ylabel('f(x)'); title('Shifted fit');

p = polyfit(xi,y,9); xi = [-5:0.1:5]; yi = polyval(p,xi);

Equivalent to interp1(x,Y,xi,’spline’)

-5 0 2 4 6 8

0

5

10

15

20

x

y

Cubic Spline

-5 0 2 4 6 8

0

5

10

15

20

x

y

7th order polynomial

• Cubic spline interpolation results in reduced fluctuation / overshoot between data points, compared to high-order polynomial fitting

Newton- Cotes

Differentiation

Interpolation

Finite difference methods

Finite Difference Methods, Interpolation, and Integration

Spline interpolation

Newton’s polynomials

Errors in FD methods

Integration

Simpson’s rules

x

f (x)

a (^) b

• Split the area bounded by the curve f (x) and the x -axis into n vertical segments (strips / rectangles)
• Calculate the area of each rectangle (height x width) and then sum the area of all rectangles to find the total area under f (x) between a and b.
• The height of each rectangle equals the value of the function at the mid- point of the rectangle
• The width of each rectangle is constant, equal to:

n

h ba

1 2 3 … (^) n

h

• How can we modify this approach to generate a better approximation for the integral?

Numerical Integration of f ( x ):

Numerical Integration / Quadrature

• Using rectangular strips is the simplest form of a more general strategy:  Replace the integrand f (x) with an approximate function whose area is easier to calculate:

     

b

a

n

b

a

I f x f x f (^) n  x  a 0 a 1 xa 2 x^2 anx^ n

 Replace the original integrand with an nth^ order polynomial function:

x

f (x)

a (^) b Zero-order polynomial (constant)

x

f (x)

a (^) b 1st-order polynomial (linear)

x

f (x)

a (^) b 2nd-order polynomial (quadratic)

The Trapezoid Rule

  ^    

b

a

b

a

I f x f 1 x f 1  x a 0 a 1 x

• Replaces the exact integrand f (x) with a 1st order polynomial function f 1 (x) :

Recall Lecture 9:

x 0 x x 1

f (x)

f (x 1 )

f 1 (x)

f (x 0 )

1 0

1 0 (^1 0) x x x x

f x f x f x f x  

  

The 1st^ order polynomial that passes through [x 0 , f (x 0 )] and [x 1 , f (x 1 )] is:

x

f (x)

a (^) b 1st-order (linear)

     ^    x a

b a

f x f a f b f a  

   1

We replace the true f ( x ) (black curve) with a 1st^ order polynomial f 1 ( x ) (red line) as an approximation for the purpose of integration

x

f (x)

a b 1st-order (linear)

      

b

a

b

a

I f x f 1 x

 x a dx

b a

f b f a I f a

b

a

 

 (^)  

  

     ^    x a

b a

f b f a f x f a  

 1  

   ^  

2

I  ba f a  f^ b

• Multiply out the second term in the square brackets, collect terms with x - dependence, integrate, rearrange:

 The trapezoid rule:

f (a)

f (b)

Example

• Numerically integrate f  x  0. 2  25 x 200 x^2  675 x^3  900 x^4  400 x^5

from a = 0 to b = 0.

(^00) 0.2 0.4 0.6 0.8 1

1

2

3

4

x

y = f(x)

(^00) 0.2 0.4 0.6 0.

1

2

3

4

x

y = f(x)

1.640533 (exact)

… a true relative error of 89.5%

2

1. 80.^20.^232 2

I  h f a  f b   

(^00) 0.2 0.4 0.6 0.8 1

1

2

3

4

x

y = f(x)

Using the Trapezoid rule gives

Trapezoid rule:

 3  2

Oh f a f b I h   The error associated with using the Trapezoid rule is O(h^3 ): 

(^00) 0.2 0.4 0.6 0.8 1

1

2

3

4

x

y = f(x) h = 0.

(^00) 0.2 0.4 0.6 0.8 1

1

2

3

4

x

y = f(x)

h = 0.

We can break down the interval from a to b into multiple smaller sections, each with smaller h, applying the Trapezoid rule to each section and summing the individual areas:

(^00) 0.2 0.4 0.6 0.8 1

1

2

3

4

x

y = f(x)^ h^ = 0.^

 This is the “multiple application” or Composite Trapezoid Rule

• The composite trapezoid rule is fine for evaluating well-behaved functions to typical levels of accuracy
• If higher accuracy is required, or if the underlying function is time consuming to evaluate, the composite trapezoid rule can become inefficient

 Can we develop more efficient numerical integration algorithms with errors that reduce more quickly as the number of segments (n) increases?

  ^    

b

a

n

b

a

I f x f x f (^) n  x  a 0 a 1 xa 2 x^2 anx^ n Use higher order polynomial functions to approximate f (x):

x

f (x)

a b

1st-order (linear)

Trapezoid rule

x

f (x)

a (^) b

2nd-order (quadratic)

Simpson’s 1/3 rule

x

f (x)

a (^) b

3rd-order (cubic)

Simpson’s 3/8 rule

Simpson’s 1/3 rule

      

b

a

b

a

I f x f 2 x f 2  x a 0 a 1 xa 2 x^2

1. Use a 2nd order polynomial function to approximate f (x):

 ^0 ^4 ^1 ^ ^2  3

f x f x f x h I    2

b a h

This can be integrated and simplified to: (^) where 

“1/3 rule”

 If we know the value of f (x) at x = a, b, and midway between a and b, ( data points) we can use Simpson’s 1/3 rule

x

f (x)

a = x 0 b = x 2

2. Approximate the integral by determining the area under this polynomial:

2 2

0 0

2 0 1 0 2 0 1

x x

x x

I  (^)  f x  (^) b  b x  x  b x  x x x dx

x 1