Math 449 - Homework 9 - Solutions
                   
                   
General note: the integral of f(x) over the interval [a,b] will be

denoted I_[a,b]f(x)dx.

                   

1) Exercise 4, section 7.1, page 362.

As we did in class for the Simpson rule, it is useful to

apply the coordinate change x = x_0 + ht, t between 0 and 1.

Then P_1(x(t)) = -f_0 (t-1) + f_1 t. Integrating from 0 to 1 and

using that dx = h dt, we obtain 

I_[x_0, x_1] P_1(x) dx = I_[0 1] (-f_0 (t-1) + f_1 t) h dt 

                       = h (f_0 + f_1)/2.

************************************************************************** 

2) Exercise 5, section 7.1, page 363.

Recall that the degree of precision of a quadrature formula is the 

largest value of n such that the truncation error E[P] is zero for

each polynomial P(x) of degree less than or equal to n. Since 

the expression E[P] is linear in P, the degree of precision is equal to

n so that E[x^i] = 0 for i < n+1 and E[x^i] is different from 0 for i=n+1.


Claim: the degree of precision for the trapezoidal rule is 1.

It is clear that this number is at least 1, since linear polynomials

agree with their first degree Lagrange interpolation approximation. Thus

all we only need to show that E[x^2] is not zero. Notice that the first

degree Lagrange polynomial with nodes x_0=0 and x_1 = 1 interpolating 

f(x) = x is x itself. Therefore,

E[x^2] = I_[0,1] x^2 dx - I_[0,1] P_1(x) dx

       = 1/3 - I_[0,1] x dx
       
       = 1/3 - 1/2
       
       =-1/6
       
which is non-zero, proving the claim.       

************************************************************************** 

3) Exercise 6, section 7.1, page 363.

Claim: the degree of precision of Simpson's rule is 3. (The Simpson's 

3/8 rule also has degree of precision = 3.)


This is just as in the previous exercise. It is clear that the degree of

precision in this case is at least 2 since we use polynomials of degree 

2 to approximate the function to be integrated. Thus it suffices to

show that E[x^3]=0 and E[x^4] is different from 0. The Lagrange polynomial

interpolating x^3 over the nodes x_0=0, x_1=1, x_2=2 is 

P_2(x) = -x(x-2) + 4x(x-1) = 3x^2 + 2x - 4. So

E[x^3] = I_[0,2] x^3 dx - I_[0,2] P_2(x) dx

       = 4 - I_[0,2](3x^2 + 2x - 4)dx
       
       = 4 - (8 + 4 - 8)
       
       = 0.
       
The second degree Lagrange polynomial interpolating x^4 over the same 

nodes is 

P_2(x) = -x(x-2) + 8x(x-1) = 7x^2 + 2x - 8. So

E[x^4] = I_[0,2] x^4 dx - I_[0,2] P_2(x) dx

       = 32/5 - I_[0,2](7x^2 + 2x - 8)dx
       
       = 32/5 - (56/3 + 4 - 16) 
       
       =-4/15,
       
which is non-zero, proving the claim.


************************************************************************** 

4) Exercise 2 (b), section 7.2, page 374.

We wish to approximate the arc-length of the curve y=sin(x) for over the

interval [0, pi/4] of x, that is, we wish to find the integral

J = I_[0, pi/4] [1 + (cos(x))^2]^(1/2) dx. Let f(x) = [1 + (cos(x))^2]^(1/2).


(i) We first do it using the composite trapezoidal rule with M=10.

h = pi/40, x_k = k*pi/40, k=0, 1, ..., 10.

T(f,h) = h[(f(0) + f(pi/4))/2 + f(x_1) + f(x_2) + ... + f(x_9)]

       = h*[(sqrt(2) + sqrt(3/2))/2 + sum(sqrt(1+cos((1:9)*h).^2))]
       
       = 1.0579

(ii) We now apply the composite Simpson rule with M=5.

h = pi/40, x_k = k*pi/40, k=0, 1, ..., 10.

S(f,h) = (h/3)[f(0) + f(pi/4)] + (2h/3)[f(x_2) + f(x_4) + ... + f(x_8)] + 
   
                               + (4h/3)[f(x_1) + f(x_3) + ... + f(x_9)]
                               
       = (h/3)*[sqrt(2) + sqrt(3/2) + 2*[sum(sqrt(1+cos((2:2:8)*h).^2))] + 
       
                                + 4*[sum(sqrt(1+cos((1:2:9)*h).^2))] ] 
                                
       = 1.0581
       
(The correct value to 5 decimal places is 1.05809.)                                                      


************************************************************************** 

5) Exercise 8 (b), section 7.2, page 375.

Let f(x) = 1/(5-x) over the interval [2, 3] and write, for an M to be

determined, h = 1/M, x_k = 2 + kh, k = 0, 1, ..., M.

I_[2,3]f(x) dx = T(f,h) + E(f,h)

were 

T(f,h) = (h/2)(f(2) + f(3)) + h[f(x_1) + ... + f(x_(M-1))]

E(f,h) = -f''(c)h^2/12

for some c in [2,3]. Thus |f''(c)| is less than or equal to 1/4 and

|E(f,h)| is less than or equal to  1/(48 M^2).

For an accuracy of 5*10^(-9), we set

1/(48 M^2) < 5*10^(-9)

so that M = 2042 will do. The corresponding value of h is 0.000489716.

************************************************************************** 

6) Algorithms and Programs 2 (b), section 7.2, page 377.

The value of the integral in that problem is bounded from below by pi/4.

Thus the relative error can be estimated by 

|E(f,h)|/(pi/4) < (pi/4)|f^(4)(c)|h^4/(180 pi/4) = |f^(4)(c)|h^4/180.

A little thought should convince you that the fourth derivative of f(x)

cannot be greater than 1800 in absolute value. (This is rather crude,

but I want to avoid the work of computing 4 derivatives of f(x).)

So the relative error is crudely estimated as less than 

10*h^4 =10*(pi/4*M)^4. Setting this number less than 5*10^(-11) gives

M = 526 is big enough. Using the program below for this M gives

the Simpson rule approximation 1.058095501392.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s=simprl(f,a,b,M)
%Input   - f is the integrand input as a string 'f'
%        - a and b are the lower and upper limits of integration
%        - M is the number of subintervals
%Output  - s is the Simpson rule sum
h=(b-a)/(2*M);
s1=0;
s2=0;
for k=1:M
    x=a+h*(2*k-1);
    s1=s1+feval(f,x);
end
for k=1:(M-1)
    x=a+h*2*k;
    s2=s2+feval(f,x);
end
s=h*(feval(f,a) + feval(f,b) + 4*s1 + 2*s2)/3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

**************************************************************************