Ma 449: Numerical Applied Mathematics
		  Model Solutions to Homework Assignment 8
			     Prof. Wickerhauser



Exercise 4 of Section 7.1, p.362.

   The integral of f0*(x-x1)/(x0-x1)+f1*(x-x0)/(x1-x0) is

	   f0*(x-x1)^2/ 2(x0-x1)+f1*(x-x0)^2/2(x1-x0)    on [x0,x1]

	   = [f1*(x-x0)^2 - f0*(x-x1)^2] / 2(x1-x0)    on [x0,x1]

	   = [f1*(x1-x0)^2 + f0*(x0-x1)^2] / 2(x1-x0)

	   = [f1*(x1-x0)^2 + f0*(x1-x0)^2] / 2(x1-x0)

	   = ( [f1+f0]/2 ) * (x1-x0),

   which is the trapezoid rule.




Exercise 5 of Section 7.1, p.363.

   Write I for the exact value of the integral and T for the trapezoid rule
   approximation.  Then, for various polynomials p=p(x):

   p0(x)=1:   I(p0) = 1*(x1-x0)
              T(p0) = ([1+1]/2) *(x1-x0)              ==> T(p0)=I(p0)

   p1(x)=x:   I(p1) = [x1^2-x0^2]/2 = ([x1+x0]/2)*(x1-x0)
	      T(p1) = ([x0+x1]/2) *(x1-x0)                ==> T(p1)=I(p1)

   p2(x)=x^2:  I(p2) = [x1^3-x0^3]/3 = ([x1^2+x0*x1+x0^2]/3)*(x1-x0)
	       T(p2) = ([x0^2+x1^2]/2) *(x1-x0)              ==> T(p2)!=I(p2)

   Hence the degree of precision of the trapezoid rule is 1.



Exercise 6 of Section 7.1, p.363.

Write I for the exact value of the integral and S for the Simpson rule
approximation.  Then, for various polynomials p=p(x):

p0(x)=1:   I(p0) = 1*(x1-x0)
	   S(p0) = ([1+4*1+1]/6) *(x1-x0)        ==> S(p0)=I(p0)

p1(x)=x:   I(p1) = [x1^2-x0^2]/2 = ([x1+x0]/2)*(x1-x0)
	   S(p1) = ([x0+4*(x0+x1)/2+x1]/6) *(x1-x0)
                 = ([x0+x1]/2) *(x1-x0)          ==> S(p1)=I(p1)

p2(x)=x^2:  I(p2) = [x1^3-x0^3]/3 = ([x1^2+x0*x1+x0^2]/3)*(x1-x0)
	    S(p2) = ([x0^2+4*(x0+x1)^2/4+x1^2]/6) *(x1-x0)
                  = ([2x0^2+2x0*x1+2x1^2]/6) *(x1-x0)
                  = ([x0^2+x0*x1+x1^2]/3) *(x1-x0)    ==> S(p2)=I(p2)

p3(x)=x^3:  I(p3) = [x1^4-x0^4]/4 = ([x1^3+x0*x1^2+x0^2*x1+x0^3]/4)*(x1-x0)
	    S(p3) = ([x0^3+4*(x0+x1)^3/8+x1^3]/6) *(x1-x0)
                  = ([x0^3+x1*x0^2+x1^2*x0+x1^3]/4)*(x1-x0) ==> S(p3)=I(p3)

p4(x)=x^4:  I(p4) = [x1^5-x0^5]/5
                  = ([x1^4+x0*x1^3+x0^2*x1^2+x0^3*x1+x0^4]/5)*(x1-x0)
	    S(p4) = ([x0^4+4*(x0+x1)^4/16+x1^4]/6) *(x1-x0)
                  = ([(5/4)x0^4 + lower-order x0 terms]/6)*(x1-x0)
				 ==> S(p4)!=I(p4)

Hence the degree of precision of Simpson's rule is 3.



Exercise 2(b) of Section 7.2, p.374.

   C PROGRAM:
   /* Exercise 2(b) i and ii, chapter 7.2, p.374 */
   #include <math.h>
   #define PI 3.14159265358979323
   typedef double real;
   typedef real (*integrand)(real); 
   /* Compute the length element for the part (b) function sin(x) */
   real
   length_element(real x) {
     real c;
     c = cos(x);			/* derivative of sin(x) */
     return sqrt(1+c*c);
   }
   /* Formula from Th 7.2, Eq 1a, p.364: */
   real
   composite_trapezoid_rule( integrand f, real a, real b, int m )
   {
     real h, sum, x;     int k;
     h = (b-a)/m;
     sum = (f(a)+f(b))/2;
     for(k=1; k<m; k++) {    x = a+k*h;    sum += f(x);  }
     sum *= h;
     return sum;
   }
   /* Formula from Th 7.3, Eq 4a, p.366: */
   real
   composite_Simpson_rule( integrand f, real a, real b, int m )
   {
     real h, sum, x;
     int k;
     h = (b-a)/(2*m);
     sum = (f(a)+f(b));
     for(k=1; k<m; k++) {    x = a+2*k*h;    sum += 2*f(x);  }
     for(k=0; k<m; k++) {    x = a+(2*k+1)*h;    sum += 4*f(x);  }
     sum *= h/3;
     return sum;
   }
   int
   main(void)
   {
     printf("Ex 7.2-2(b):  length of f(x)=sin(x) over [0,pi/4]\n");
     printf("Composite trapezoid rule value, M=10: %20.14f\n",
	    composite_trapezoid_rule(length_element,0,PI/4,10) );
     printf("Composite Simpson rule value, M=5: %20.14f\n",
	    composite_Simpson_rule(length_element,0,PI/4,5) );
     return 0;
   }

   OUTPUT:
   % ./a.out
   Ex 7.2-2(b):  length of f(x)=sin(x) over [0,pi/4]
   Composite trapezoid rule value, M=10:     1.05788556554646
   Composite Simpson rule value, M=5:     1.05809581806702


Exercise 8(b) of Section 7.2, p.375.

   The error is bounded by e=M2*(b-a)*h^2/12, where M2 is the maximum value of
   |f''(c)| on [a,b], and h=(b-a)/m.  But for f(x)=1/(5-x) on [2,3], we have
   f''(x) = 2/(5-x)^3, which is no bigger than 1/4 on [2,3].  Also, (b-a)=1
   and h=1/m for this interval.  Thus, to satisfy e<5e-9, we must have

	   h < sqrt( 12*e / M2 ) = sqrt( 240 e -9 ) ~ 5 e -4,
	   m > sqrt( M2 / 12*e ) ~ 2 e 3,

   so m=2200 should work.



Algorithm 2(b) of Section 7.2, p.377.

   The composite Simpson rule error is bounded by e=M4*(b-a)*h^4/180, where
   M4 is the maximum value of |f''''(c)| on [a,b], and h=(b-a)/2m.  But for
   f(x)=sqrt(1+cos(x)^2) on [0,pi/4], we have

    f''''(x) = (  7290 + 11986 Cos[2 x] + 1128 Cos[4 x] + 45 Cos[6 x] + 

					       9/2               9/2
		  30 Cos[8 x] + Cos[10 x]) / (2    (3 + Cos[2 x])   ),


   which is no bigger than 2.2 on [0,pi/4].  [This can be estimated with
   Matlab's symbolic algebra and plotting routines.]

   Also, (b-a)=pi/4 and h=pi/8m for this interval.  Thus, to satisfy e<5e-12
   for 11 decimal places, we must have

	   h^4 <  920*e / M4*pi ==> h < .005,
	   m > pi/8h ~ 77

   so m=80 should work.

   C PROGRAM:
   /* Algorithm 2 for Exercise 2(b) ii, chapter 7.2, p.377 */
   #include <math.h>
   typedef double real;
   typedef real (*integrand)(real); 
   #define PI 3.14159265358979323
   /* Compute the length element for the part (b) function sin(x) */
   real
   length_element(real x) {
     real c;
     c = cos(x);			/* derivative of sin(x) */
     return sqrt(1+c*c);
   }
   /* Formula from Th 7.3, Eq 4a, p.376: */
   real
   composite_Simpson_rule( integrand f, real a, real b, int m )
   {
     real h, sum, x;  int k;
     h = (b-a)/(2*m);
     sum = (f(a)+f(b));
     for(k=1; k<m; k++) {    x = a+2*k*h;    sum += 2*f(x);  }
     for(k=0; k<m; k++) {    x = a+(2*k+1)*h;    sum += 4*f(x);  }
     sum *= h/3;
     return sum;
   }
   int
   main(void)
   {
     printf("Alg 7.2-2 for Ex 2(b):  length of f(x)=sin(x) over [0,pi/4]\n");
     printf("Composite Simpson rule value, M=80: %20.14f\n",
	    composite_Simpson_rule(length_element,0,PI/4,80) );
     return 0;
   }

   OUTPUT:
   % ./a.out            
   Alg 7.2-2 for Ex 2(b):  length of f(x)=sin(x) over [0,pi/4]
   Composite Simpson rule value, M=80:     1.05809550139739