4.1    performance of Trapezoid Rule for concave/convex Functions

(a).  prove that if a function is concave downward, then the trapezoid rule underesti- mates the integral, i.e., it gives a value smaller than the exact solution.

(b).   similarly,  prove  that  if a  function is convex upward,  then  the trapezoid rule overestimates the integral.

4.2    A General property

show that there exists coe伍cients (山0, 山1 , . . . , 山n) depending on (①0, ① 1 , . . .  , ①n) and (a, b) such that

labp(①) d① =  山ip(①i)

for all polynomials p of degree 三  n.

( Hint: use the Lagrange form of the interpolating polynomial.)

4.3    Trapezoid and simpson,s Methods

Given function f (①) = e-① , we study diferent numerical approximations to the integral

l00(0.)8 f (①) d① .

we will use the values of f (①) at the points O.O, O.2, O.4, O.6, O.8.  Generate the data set before you start the numerical integration. use 6-digits accuracy.

(a). write out the trapezoid rule and compute the numerical integration with 6 digits.

(b). write out the simpson,s rule and compute the numerical integration with 6 digits.

(c). what is the exact value of the integral？  what is the absolute error by using trapezoid and simpson,s rule？which method is better？

(d).  The error formula for the trapezoid rule with n + 1 points yields

b — a    2      地地                       b — a

12                             n   ,

for some ξ e (a, b). The error for simpson,s rule with (2n+1) points yields

Es(f ; h) = — h4 f(4)(ξ地 ),    h =  ,

for some ξ地  = (a, b).  If we wish the absolute value of the error to be smaller than 1O-4 , how many points would be needed for each method？

4.4    simpson,s rule.

consider the simpson,s rule for the integral    l01 f (①) d①,   with the given data set:

The plots of the derivatives of f are given below:

1st derivative

1.0

0.5

0.0

0.0    0.2    0.4    0.6    0.8    1.0

2nd derivative

2.0

1.8

1.6

1.4

1.2

1.0

0.8

0.6

0.4

0.2

0.0    0.2    0.4    0.6    0.8    1.0

3rd derivative

- 0.0

- 0.5

- 1.0

- 1.5

- 2.0

- 2.5

- 3.0

- 3.5

0.0    0.2    0.4    0.6    0.8    1.0

4th derivative

- 1.0

- 1.5

- 2.0

- 2.5

- 3.0

- 3.5

- 4.0

0.0    0.2    0.4    0.6    0.8    1.0

(a). compute the numerical integration using simpson,s rule with h = O.25.

(b) If the error tolerance is 1O-6 , how many points must one use for simpson,s rule with uniform grid？

4.5    Trapezoid Rule in Matlab

preparation:   use the help in Matlab to learn how to use the function feval.

write a Matlab function which computes the integral by trapezoid rule. your function should be used by the following command in Matlab command window:

>> v=trape处oid（,funⅠtg,,a,b,n）

where funⅠtg.m is the name of the ile of the function f (①), and a,b is the interval, and n is the number of sub-intervals (i.e., n+1 will be the number of points.)

Test your function on computing the integral

l00(0.)8 e-① d① .

write a script that computes for n = 4, 8, 16, 32, 64, 128.  compute also the absolute error for each n,  and make a plot of the absolute error against  n.   (use  loglog  to plot.) How does the error change when n is doubled？ Do you expect this from the error estimate？ write your comments.

what to hand in？   submit your function in the ile trape处oid.m, your script ile, the plot of error, and your comments.

4.6    simpson,s Rule in Matlab

write a Matlab function which computes the integral by simpson,s rule. your function should be used by the following command in Matlab command window:

>> v=simpson（,funⅠtg,,a,b,n）

where funⅠtg.m is the name of the ile of the function f (①), and a,  b is the interval, and n is the number of sub-intervals (i.e., 2n+1 will be the number of points.)

Test your function on problem 1. write ascript that computes, for n = 2, 4, 8, 16, 32, 128, the absolute error for each n, and a plot (with loglog) of the absolute error against n. How does the error change when n doubled？ compare the results with trapezoid rule and comment.

what to hand in？   your function in simpson.m, your script ile, the plot of error,

and your comments.

4.7    simpson,s Rule over Non-uniform Mesh

For a given grid a = ①0   < ① 1   <  . . .  <  ①N-1   < ①N  = b, not necessarily uniform, the simpson,s rule for the numerical integration is

. [f (①i) + 4f (①i+

) + f (①i+1)] .

(1)

Here ① = ①i + ①i+1) is the mid-point of the interval [①i, ①i+1].

(a). Implement the simpson,s rule (1) on an arbitrary mesh:

function  Ⅰ  =  simpsonG（f,  xg）

Description of the parameters:

— f : function handle or string

— xg: row vector (①0, ① 1 , . . .  , ①N-1, ①N).

Note that you do not have to input a and b, since ①0  = a and ①N  = b.

(b).  Let N be mesh parameters.  we construct a mesh over the interval [o,1] by setting

从k = 从k-1 + 2-k ,    k = 1, 2,    , N — 1

write a Matlab function “GenMesh” that generates this mesh:

function  xg  = GenMesh（N）

The output should be a row vector representing the mesh.

(c).  Test the code on the function f (从) = e① — 从, integrated over [o, 1], for N = 4, 8, 16. show the answers from your Matlab functions.

The exact value of the integral can be easily worked out. please work it out. compute the error of the simpson,s method, for N = 4, 8, 16.

submit the script ile and the output of the code.

4.8    Trapezoid rule and Romberg algorithm.

(a) compute the trapezoid rule approximations for the integral  ' 11 3从2 d从 with n = 1 and n = 2.

(b) Apply Romberg algorithm to the results in part (a) to get a better value.  Did you get the exact value？ would you like to comment on it？

4.9    Romberg Algorithm in Matlab

preparation:   use helpdesk in Matlab to learn how to use the functions feval and

quad.

(a). write a Matlab function that computes the Romberg integration.  one should be able to call the function by:

>> R=romberg（,f,,a,b,n）

where f is the name of the function where f (从) is implemented, and a and bdeines the integrating interval, and n is the size of your Romberg table.  The function should return the whole Romberg table.  The best approximation of the interval would be the value in R（n,n）.

You may follow the pseudo-code in the lecture notes, or implement your own.

You may check your code against the result of the simulation in section 4.8, to make sure that your code works.

(b).  use your romberg to compute the integrals

i)

ii)

l0 π sin(①) d①

l01 ^① d①

(= 2)

(= 2/3)

compute also the errors. The exact values of the integrals are given in the paren- theses above.  print the errors along the diagonal of the table, and note how it changes along the diagonal of the table.  use format  short  e in Matlab to display the error data.

(c). Explain why Romberg algorithm works poorly for the integral in ii).

(d). use Matlab functions quad and quadl to compute the integrals in b).  use  1e-9 as tolerance for both integrations. Mark your observation.

what to hand in？   The Matlab ile romberg.m, the ile for your function f.m, a script ile that does b) and d), and the Romberg tables your get in b).

4.10    computing π with Romberg

The transcendental number π is the area of a circle with radius 1.  show that

π = 8 l01/^2  (^1 - ①2  - ①) d① .

You may use a diagram to see this more clearly.

computing the value of π using your Romberg algorithm, with a Romberg triangle of size at least 6.

submit your Matlab function/script, and the output result.

4.11    Numerical Integration and Extrapolation

consider the function

! sint

t  O,

t = O.

Note that f is a continuous function for all t. Let

J = l01 f (t)dt.

(a).  compute approximations to J by trapezoid rule with 1, 3 and 9 equal intervals. Try to use fewest possible time of computing f (t), similar to what we did with recursive trapezoid rule. show which formula you use.

(b).  Derive a Richardson extrapolation algorithm which uses trapezoid approximations from a), to obtain better approximations to  J.   show  which  formula  you  use. Generate the corresponding Romberg triangle.

4.12    Gaussian Quadrature and Beyond

(a).    consider the Gaussian Quadrature rule with 4 points on the interval [-1, 1],

l 11 f (①) d① 必 a1 f (① 1) + a2 f (①2) + a3 f (①3) + a4 f (①4)

where

① 1     =   -  (3 - 4^0.3),

①2     =   - (3 + 4^0.3),

①3     =    (3 - 4^0.3),

①4     =    (3 + 4^0.3),

where

a1     =    +   ,

a2     =    -   ,

a3     =    +   ,

a4     =    -  .

show that the rule is exact for all polynomials of degree 参 7.

(b).    construct a rule of the form

l 11 f (①) d① 必 a1 f (-0.5) + a2 f (0) + a3 f (0.5)

that is exact for all polynomials of degree 参 2；that is, determine the values for a1 , a2 , a3 .

(c).    Find positive values 山 and a such that the following rule is exact for polynomials of degree 三 3

l 11 f (①) d① 必 山f (-a) + 山f (a).

4.13    More Fun practice

(a). Determine a formula of the form

loh f (①) d① 必 山of (o) + 山1f (h) + 山2f\\ (o) + 山3f\\ (h)

that is exact for polynomials of as high a degree as possible.

(b). Derive a numerical integration formula of the form

laah(+)h f (①) d① 必 Af (a) + Bf\ (a - h) + cf\\ (a + h)

that is exact for polynomials of as high a degree as possible.

4.14    Error Formula for simpson,s Rule and simpson,s

(a). consider the simpson,s rule, which can be written as

Rule； (Bonus)