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



Exercise 4 of Section 7.3, p.390.

   First we note that the integral of t^k from 0 to 4 is
   4^(k+1)/(k+1) and takes the following values:

   k	integral of t^k from 0 to 4
   -	---------------------------
   0	4
   1	8
   2	64/3
   3	64
   4	1024/5

   Next, set the values of these integrals equal to a linear
   combination of function values 0^k, 1^k, 2^k, 3^k, and 4^k,
   with unknown weights w0,w1,w2,w3,w4, and solve for the w's:

   MATLAB CODE:
    <>t=[0 1 2 3 4];
    <>a=[ones(t);t;t.*t;t.*t.*t;t.*t.*t.*t];
    <>b=[4;8;64/3;64;1024/5];
    <>long;
    <>w=a\b

    W     =

       .311111111111114
      1.422222222222215
       .533333333333342
      1.422222222222218
       .311111111111112

    <>(45/2)*w

    ANS   =

      7.000000000000060
     31.999999999999833
     12.000000000000185
     31.999999999999897
      7.000000000000024

   Except for the round-off error, we recognize that
   W = (2/45)*[ 7 32 12 32 7 ] 






Algorithm 3 of Section 7.3, p.392.

   We do this by Romberg integration:
   C PROGRAM:
   /* Section 7.3, Algorithm 3, done with Romberg integration */
   #include <assert.h>
   #include <math.h>
   #include <stdio.h>
   #define PI 3.14159265358979323
   #define JMAX 15			/* Recursion depth limit. */
   /* Function to integrate to get Phi: */
   double f(double t) {  return exp(-t*t); }
   /* Recursive composite trapezoid rule increment on 2^J segments of [a,b]: */
   double rcti(double a, double b, int j)
   {
     double h, sum;  int n, k;
     assert(a<b);  assert(j>=0);  assert(j<30);
     n = 1<<j;  h = (b-a)/n;
     if(j==0) sum = (f(a)+f(b))/2;
     else { sum = 0; for(k=1; k<n; k+=2) sum += f(a+k*h);  }
     return h*sum;
   }
   main(void)
   {
     double r[JMAX][JMAX], phi;
     int i, j, k;
     double a=0.0, b[8]={0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0};
     double err, tol=5.0e-14;

     for(i=0; i<8; i++) {
       r[0][0] = rcti(a,b[i],0);	/* Initial trapezoid rule value */
       for(j=1; j<JMAX; j++) {
	 r[j][0]=0.5*r[j-1][0]+rcti(a,b[i],j); /* recursive trapezoid rule */
	 for(k=1; k<=j; k++)		/* Romberg integration */
	   r[j][k] = r[j][k-1] + (r[j][k-1] - r[j-1][k-1])/((1<<(2*k))-1.0);
	 err = fabs( r[j][j] - r[j-1][j-1] );
	 phi = 0.5 + r[j][j]/sqrt(2*PI); /* adjustment, to get Phi */
	 printf("phi[%3.1f]=%20.14f    err=%8.2e, J=%d\n", 
		b[i], phi, err, j );
	 if(err<tol) break;	/* escape j-loop, go to next value b[i] */
       }
     }
     return 0;
   }

   OUTPUT:
   phi[0.5]=    0.69147333953143    err=9.33e-03, J=1
   phi[0.5]=    0.69146244243255    err=2.73e-05, J=2
   phi[0.5]=    0.69146246128594    err=4.73e-08, J=3
   phi[0.5]=    0.69146246127401    err=2.99e-11, J=4
   phi[0.5]=    0.69146246127401    err=6.05e-15, J=5

   phi[1.0]=    0.84152905199630    err=5.28e-02, J=1
   phi[1.0]=    0.84134391691401    err=4.64e-04, J=2
   phi[1.0]=    0.84134474699460    err=2.08e-06, J=3
   phi[1.0]=    0.84134474606877    err=2.32e-09, J=4
   phi[1.0]=    0.84134474606854    err=5.82e-13, J=5
   phi[1.0]=    0.84134474606854    err=3.33e-16, J=6

   phi[1.5]=    0.93325240117164    err=9.25e-02, J=1
   phi[1.5]=    0.93320473311640    err=1.19e-04, J=2
   phi[1.5]=    0.93319266781615    err=3.02e-05, J=3
   phi[1.5]=    0.93319279907011    err=3.29e-07, J=4
   phi[1.5]=    0.93319279873091    err=8.50e-10, J=5
   phi[1.5]=    0.93319279873114    err=5.93e-13, J=6
   phi[1.5]=    0.93319279873114    err=8.88e-16, J=7

   phi[2.0]=    0.97360538166373    err=5.18e-02, J=1
   phi[2.0]=    0.97744077645329    err=9.61e-03, J=2
   phi[2.0]=    0.97724742829474    err=4.85e-04, J=3
   phi[2.0]=    0.97724987667571    err=6.14e-06, J=4
   phi[2.0]=    0.97724986804334    err=2.16e-08, J=5
   phi[2.0]=    0.97724986805182    err=2.13e-11, J=6
   phi[2.0]=    0.97724986805182    err=5.33e-15, J=7
   phi[2.5]=    0.97794455102129    err=1.07e-01, J=1
   phi[2.5]=    0.99470541419443    err=4.20e-02, J=2
   phi[2.5]=    0.99377636765482    err=2.33e-03, J=3
   phi[2.5]=    0.99379039175009    err=3.52e-05, J=4
   phi[2.5]=    0.99379033461630    err=1.43e-07, J=5
   phi[2.5]=    0.99379033467424    err=1.45e-10, J=6
   phi[2.5]=    0.99379033467422    err=2.82e-14, J=7

   phi[3.0]=    0.96072225573847    err=3.62e-01, J=1
   phi[3.0]=    1.00099655791265    err=1.01e-01, J=2
   phi[3.0]=    0.99861513281805    err=5.97e-03, J=3
   phi[3.0]=    0.99865019272069    err=8.79e-05, J=4
   phi[3.0]=    0.99865010209893    err=2.27e-07, J=5
   phi[3.0]=    0.99865010196812    err=3.28e-10, J=6
   phi[3.0]=    0.99865010196837    err=6.33e-13, J=7
   phi[3.0]=    0.99865010196837    err=4.44e-16, J=8

   phi[3.5]=    0.93453913906814    err=6.65e-01, J=1
   phi[3.5]=    1.00349313493132    err=1.73e-01, J=2
   phi[3.5]=    0.99974222565594    err=9.40e-03, J=3
   phi[3.5]=    0.99976690382516    err=6.19e-05, J=4
   phi[3.5]=    0.99976737360765    err=1.18e-06, J=5
   phi[3.5]=    0.99976737091816    err=6.74e-09, J=6
   phi[3.5]=    0.99976737092097    err=7.03e-12, J=7
   phi[3.5]=    0.99976737092096    err=2.00e-15, J=8

   phi[4.0]=    0.91002665111997    err=9.73e-01, J=1
   phi[4.0]=    1.00339140917081    err=2.34e-01, J=2
   phi[4.0]=    1.00007534563051    err=8.31e-03, J=3
   phi[4.0]=    0.99996535049684    err=2.76e-04, J=4
   phi[4.0]=    0.99996834162712    err=7.50e-06, J=5
   phi[4.0]=    0.99996832874534    err=3.23e-08, J=6
   phi[4.0]=    0.99996832875817    err=3.22e-11, J=7
   phi[4.0]=    0.99996832875817    err=8.22e-15, J=8

   So our final values are:
   phi[0.5]=    0.69146246127401    err=6.05e-15, J=5
   phi[1.0]=    0.84134474606854    err=3.33e-16, J=6
   phi[1.5]=    0.93319279873114    err=8.88e-16, J=7
   phi[2.0]=    0.97724986805182    err=5.33e-15, J=7
   phi[2.5]=    0.99379033467422    err=2.82e-14, J=7
   phi[3.0]=    0.99865010196837    err=4.44e-16, J=8
   phi[3.5]=    0.99976737092096    err=2.00e-15, J=8
   phi[4.0]=    0.99996832875817    err=8.22e-15, J=8

   Note: Phi(x) = 0.5*(1+Erf[x/Sqrt[2]]), where Erf[] is the
	better-known "error function".  A table of values of
	Phi(x) produced by using Mathematica's built-in Erf is:

      Phi(0.5)= 0.69146246127401
      Phi(1.0)= 0.84134474606854 
      Phi(1.5)= 0.93319279873114
      Phi(2.0)= 0.97724986805182
      Phi(2.5)= 0.99379033467422
      Phi(3.0)= 0.99865010196837
      Phi(3.5)= 0.99976737092096
      Phi(4.0)= 0.99996832875816







Algorithm 1(c) of Section 7.4, p.399.

   C PROGRAM:
   /* Adaptive quadrature using Simpson's rule */
   #include <math.h>
   #include <stdio.h>
   #include <stdlib.h>
   int depth;			/* used to track recursion depth */
   double f(double x){ return 1/sqrt(x); }
   double simpson(double a, double b)
   {
     return (b-a)*(f(a)+4*f((a+b)/2)+f(b))/6;
   }
   double rsrule(double a, double b, double sum, double tol, int descend)
   {
     double sr1, sr2, c, err;

     if(depth>descend) depth=descend; /* track recursion depth  */
     if(!descend) return sum;	/* recursion limit is reached */
     c = (a+b)/2;			/* else descend if necessary */
     sr1 = simpson(a,c);  sr2 = simpson(c,b);
     err = fabs(sr1+sr2 - sum);	/* test the current error */
     if(err<tol)			/* accuracy is achieved */
       sum = sr1+sr2;
     else				/* split the current interval */
       sum =
	 rsrule(a,c,sr1,tol/2,descend-1)
	 + rsrule(c,b,sr2,tol/2,descend-1);
     return sum;
   }
   main(char argc, char **argv)
   {
     double tol0, a0=0.04, b0=1;
     double sum, h;

     tol0 = atof(argv[1]);	/* get tolerance from the command line */
     sum = simpson(a0, b0);	/* start the recursion */
     depth = 20;		/* allow a reasonably deep maximum level */
     h = (b0-a0)*pow(0.5,(double)depth); /* smallest allowed subinterval */
     sum = rsrule(a0,b0,sum,tol0,depth); /* recursive Simpson rule */
     h *= pow(2.0,(double)depth);	/* smallest interval actually used */
     printf("sum=%.16f [tol=%.0e]  smallest h=%.2e (%d from bottom)\n",
	    sum, tol0, h, depth);
     return 0;
   }

   OUTPUT:
   sum=1.6000000000000001 [tol=1e-15]  smallest h=7.32e-06 (3 from bottom)
   sum=1.6000000000000003 [tol=1e-14]  smallest h=2.93e-05 (5 from bottom)
   sum=1.6000000000000023 [tol=1e-13]  smallest h=5.86e-05 (6 from bottom)
   sum=1.6000000000000241 [tol=1e-12]  smallest h=1.17e-04 (7 from bottom)
   sum=1.6000000000002295 [tol=1e-11]  smallest h=1.17e-04 (7 from bottom)
   sum=1.6000000000019940 [tol=1e-10]  smallest h=2.34e-04 (8 from bottom)
   sum=1.6000000000213106 [tol=1e-09]  smallest h=4.69e-04 (9 from bottom)
   sum=1.6000000002319235 [tol=1e-08]  smallest h=9.37e-04 (10 from bottom)
   sum=1.6000000023097760 [tol=1e-07]  smallest h=1.87e-03 (11 from bottom)
   sum=1.6000000197334954 [tol=1e-06]  smallest h=3.75e-03 (12 from bottom)
   sum=1.6000001709258127 [tol=1e-05]  smallest h=7.50e-03 (13 from bottom)
   sum=1.6000023008540003 [tol=1e-04]  smallest h=1.50e-02 (14 from bottom)
   sum=1.6000345963877050 [tol=1e-03]  smallest h=3.00e-02 (15 from bottom)
   sum=1.6000642999974843 [tol=1e-02]  smallest h=6.00e-02 (16 from bottom)
   sum=1.6030655628075938 [tol=1e-01]  smallest h=2.40e-01 (18 from bottom)






Exercise 9 of Section 7.5, p.406.

   Since both the integral and the 3-point Gauss-Legendre quadrature
   rule are zero for the odd functions x, x^3, and x^5, it is not
   necessary to compute those functions' integrals.

   For the even functions, we have:

	indefinite   definite        3-point Gauss-Legendre
   f(x)	integral     integral        quadrature rule value
   ----	----------   --------	     ---------------------
    1       x            2                     2		
   x^2	x^3 / 3	        2/3          (5*0.6+0+5*0.6)/9 = 2/3
   x^4	x^5 / 5	     2/5 = 0.4       (5*(.6)^2+0+5(.6)^2)/9 = 0.4
   x^6	x^7 / 7	     2/7 ~ 0.285714  (5*(.6)^3+0+5(.6)^3)/9 = 0.24

   Hence the degree of precision of this quadrature rule is 5.







Problem H09.1.  Write a program to compute the integral of a function using
successive refinement of h with the composite trapezoidal rule.  Let the
desired accuracy be specified by a parameter in the program.

   This program is very similar to the one developed for Algorithm 1(c)
   of Section 7.4, p.399.

   C PROGRAM:
   /* Adaptive quadrature using trapezoid rule */
   #include <math.h>
   #include <stdio.h>
   #include <stdlib.h>
   int depth;			/* used to track recursion depth */
   double f(double x){ return 1/sqrt(x); }
   double trapezoid(double a, double b)
   {
     return (b-a)*(f(a)+f(b))/2;
   }
   double rtrule(double a, double b, double sum, double tol, int descend)
   {
     double tr1, tr2, c, err;
     if(depth>descend) depth=descend; /* track recursion depth  */
     if(!descend) return sum;	/* recursion limit is reached */
     c = (a+b)/2;			/* else descend if necessary */
     tr1 = trapezoid(a,c);  tr2 = trapezoid(c,b);
     err = fabs(tr1+tr2 - sum);	/* test the current error */
     if(err<tol)			/* accuracy is achieved */
       sum = tr1+tr2;
     else				/* split the current interval */
       sum =
	 rtrule(a,c,tr1,tol/2,descend-1)
	 + rtrule(c,b,tr2,tol/2,descend-1);
     return sum;
   }
   main(char argc, char **argv)
   {
     double tol0, a0=0.04, b0=1;
     double sum, h;

     tol0 = atof(argv[1]);	/* get tolerance from the command line */
     sum = trapezoid(a0, b0);	/* start the recursion */
     depth = 20;		/* allow a reasonably deep maximum level */
     h = (b0-a0)*pow(0.5,(double)depth); /* smallest allowed subinterval */
     sum = rtrule(a0,b0,sum,tol0,depth); /* recursive trapezoid rule */
     h *= pow(2.0,(double)depth);	/* smallest interval actually used */
     printf("sum=%.16f [tol=%.0e]  smallest h=%.2e (%d from bottom)\n",
	    sum, tol0, h, depth);
     return 0;
   }

   OUTPUT:
   sum=1.6091158094120415 [tol=1e-01]  smallest h=3.00e-02 (15 from bottom)
   sum=1.6020144989836176 [tol=1e-02]  smallest h=7.50e-03 (13 from bottom)
   sum=1.6001809524884210 [tol=1e-03]  smallest h=1.87e-03 (11 from bottom)
   sum=1.6000182313760580 [tol=1e-04]  smallest h=4.69e-04 (9 from bottom)
   sum=1.6000016835774629 [tol=1e-05]  smallest h=2.34e-04 (8 from bottom)
   sum=1.6000001858215835 [tol=1e-06]  smallest h=5.86e-05 (6 from bottom)
   sum=1.6000000182655667 [tol=1e-07]  smallest h=1.46e-05 (4 from bottom)
   sum=1.6000000016960227 [tol=1e-08]  smallest h=7.32e-06 (3 from bottom)
   sum=1.6000000001856700 [tol=1e-09]  smallest h=1.83e-06 (1 from bottom)
   sum=1.6000000000202697 [tol=1e-10]  smallest h=9.16e-07 (0 from bottom)
   sum=1.6000000000053904 [tol=1e-11]  smallest h=9.16e-07 (0 from bottom)
   sum=1.6000000000043675 [tol=1e-12]  smallest h=9.16e-07 (0 from bottom)
   sum=1.6000000000043304 [tol=1e-13]  smallest h=9.16e-07 (0 from bottom)
   sum=1.6000000000043304 [tol=1e-14]  smallest h=9.16e-07 (0 from bottom)

   Asking for more accuracy clearly gains nothing.