Math 449 - Homework 10 - Solutions
                   
                   
The integral of f(x) over the interval [a,b] will be denoted I_[a,b]f(x)dx. 


1) Exercise 1 (b), section 7.3, page 389.

Construct by hand a Romberg table with three rows for the integral

I_[0, 3] sin(4x) e^(-2x) dx = 0.1997146621.


The first column is 

R(0,0) = T(0) = (3/2)(f(0) + f(3))              = -0.00199505
R(1,0) = T(1) = T(0)/2 + (3/2)f(3/2)            = -0.02186444
R(2.0) = T(2) = T(1)/2 + (3/4)(f(3/4) + f(9/4)) =  0.01611754

The second column is 

R(1,1) = (4R(1,0) - R(0,0))/3 = -0.02848757
R(2,1) = (4R(2,0) - R(1,0))/3 =  0.02877821

The third column is

R(2,2) = (16R(2,1) - R(1,1))/15 = 0.03259592

So we have the table:

-0.00199505
-0.02186444 -0.02848757
 0.01611754  0.02877821  0.03259592

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

2) Exercise 4, section 7.3, page 390.

Derive Boole's rule (M = 1, h = 1) by using the method of undetermined

coefficients: Find the constants w_0, w_1, w_2, w_3, w_4 so that 

I_[0,4] g(t) dt = w_0g(0) + w_1g(1) + w_2g(2) + w_3g(3) + w_4g(4)

is exact for the five functions g(t) = 1, t, t^2, t^3, t^4.


The values of I_[0,4] g(t) dt are:

I_[0,4] 1 dt   = 4
I_[0,4] t dt   = 8
I_[0,4] t^2 dt = 64/3
I_[0,4] t^3 dt = 64
I_[0,4] t^4 dt = 1024/5

Thus we obtain the linear system:

w_0 + w_1 +   w_2 +   w_3 +    w_4 = 4
      w_1 +  2w_2 +  3w_3 +   4w_4 = 8
      w_1 +  4w_2 +  9w_3 +  16w_4 = 64/3
      w_1 +  8w_2 + 27w_3 +  64w_4 = 64
      w_1 + 16w_2 + 81w_3 + 256w_4 = 1024/5
      
The solution to this system is:

w_0 = 14/45, w_1 = 64/45, w_2 = 8/15, w_3 = 64/45, w_4 = 14/45.


Therefore, the quadrature formula is:

I_[0,4] g(t) dt = (2/45)[7g(0) + 32g(1) + 12g(2) + 32g(3) + 7g(4)]        


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

3) Algorithms and Programs 3, section 7.3, page 392.

The normal probability density function is 

                         f(t) = (1/sqrt(2pi))exp(-t^2/2)

and the cumulative distribution is a function defined by 

                         F(t) = (1/2) + (1/sqrt(2pi)) I_[0,x]exp(-t^2/2) dt.
                         
Compute values for 

F(0.5), F(1.0), F(1.5), F(2.0), F(2.5), F(3.0), F(3.5), F(4.0)

with 8 digits of accuracy.


I used program 7.4 with parameters: n=10, a=0, b=0.5, 1.0, etc., tol=10^(-9).

The function file for f(x) (written g.m below) is 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=g(x)
y=1/2 + (1/sqrt(2*pi))*exp(-x^2/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The program is shown below. To find F(0.5), for example, I used the command:

format long
[R,quad,err,h]=romber('g',0,0.5,10,10^(-9));
quad

I obtained the following:

F(0.5) = 0.691462461
F(1.0) = 0.841344746
F(1.5) = 0.933192799
F(2.0) = 0.977249868
F(2.5) = 0.993790335
F(3.0) = 0.998650102
F(3.5) = 0.999767371
F(4.0) = 0.999968329


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [R,quad,err,h]=romber(f,a,b,n,tol)
%Input  - f is the integrand input as a string 'f'
%       - a and b are upper and lower limits of integration
%       - n is the maximum number of rows in the table
%       - tol is the tolerance
%Output - R is the Romberg table
%       - quad is the quadrature value
%       - err is the error estimate
%       - h is the smallest step size used
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
while((err>tol)&(J<n))|(J<4)
    J=J+1;
    h=h/2;
    s=0;
    for p=1:M
        x=a+h*(2*p-1);
        s=s+feval(f,x);
    end
    R(J+1,1)=R(J,1)/2 + h*s;
    M=2*M;
    for K=1:J
        R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
    end
    err=abs(R(J,J)-R(J+1,K+1));
end
quad=R(J+1,J+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

4) Algorithms and Programs 1 (c), section 7.4, page 399.

Use Program 7.6 to approximate the value of I_[0.04 1] 1/sqrt(x) dx with

starting tolerance e_0 = 0.00001.


The relevant programs are shown below.

The command (after saving a function file for g(x)=1/sqrt(x)) is 

[SRmat, quad, err]=adapt('g',0.04,1,0.00001);

which gives

quad = 1.6000023

with 

err = 3.23 10^(-6)

The exact value can be easily calculated and equals 1.6


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Z=srule(f, a0, b0, tol0)
%Input  - f is the integrand input as a string 'f'
%       - a0 and b0 are the lower and upper limits of integration
%       - tol0 is the tolerance
%Output - Z is a 1x vector [a0 b0 S S2 err tol1]
h=(b0-a0)/2;
C=zeros(1,3);
C=feval(f,[a0 (a0+b0)/2 b0]);
S=h*(C(1) + 4*C(2) + C(3))/3;
S2=S;
tol1=tol0;
err=tol0;
Z=[a0 b0 S S2 err tol1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [SRmat, quad, err]=adapt(f,a,b,tol)
%Input  - f is the integrand input as a string 'f'
%       - a and b are lower and upper limits of integration
%       - tol is the tolerance
%Output - SRmat is the table of values
%       - quad is the quadrature value
%       -err is the error estimate
%Initialize vaues
SRmat = zeros(30,6);
iterating=0;
done=1;
SRvec=zeros(1,6);
SRvec=srule(f,a,b,tol);
SRmat(1,1:6)=SRvec;
m=1;
state=iterating;

while(state==iterating)
    n=m;
    for j=n:-1:1
        p=j;
        SR0vec=SRmat(p,:);
        err=SR0vec(5);
        tol=SR0vec(6);
        if (tol<=err)
            %Bisect interval, apply Simpson's rule
            %recursively, and determine error
            state=done;
            SR1vec=SR0vec;
            SR2vec=SR0vec;
            a=SR0vec(1);
            b=SR2vec(2);
            c=(a+b)/2;
            err=SR0vec(5);
            tol=SR0vec(6);
            tol2=tol/2;
            SR1vec=srule(f,a,c,tol2);
            SR2vec=srule(f,c,b,tol2);
            err=abs(SR0vec(3)-SR1vec(3)-SR2vec(3))/10;
            
            %Accuracy test
            if (err<tol)
                SRmat(p,:)=SR0vec;
                SRmat(p,4)=SR1vec(3)+SR2vec(3);
                SRmat(p,5)=err;
            else
                SRmat(p+1:m+1,:)=SRmat(p:m,:);
                m=m+1;
                SRmat(p,:)=SR1vec;
                SRmat(p+1,:)=SR2vec;
                state=iterating;
            end
        end
    end
end
quad=sum(SRmat(:,4));
err=sum(abs(SRmat(:,5)));
SRmat=SRmat(1:m,1:6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        

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