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





Exercise 4 of Section 6.1, p.335.

   (a) Formula 10 MATLAB code and its output:

    <>num = [1 -8 0 8 -1]; dx=[-2 -1 0 1 2];
    <>x1=2.3*ones(dx) +0.1*dx,  f1=exp(x1)
    X1    =    2.1000    2.2000    2.3000    2.4000    2.5000
    F1    =    8.1662    9.0250    9.9742   11.0232   12.1825
    <>x01=2.3*ones(dx) +0.01*dx,  f01=exp(x01)
    X01   =    2.2800    2.2900    2.3000    2.3100    2.3200
    F01   =    9.7767    9.8749    9.9742   10.0744   10.1757
    <>h=0.1; den=12*h; df1=(f1*num')/den;
    <>long; df1, short; err=df1-exp(2.3)
    DF1   =    9.974149167936689		 ERR   =   -3.3287E-05
    <>h=0.01; den=12*h; df01=(f01*num')/den;
    <>long; df01, short; err=df01-exp(2.3)
    DF01  =    9.974182451489710		 ERR   =   -3.3250E-09

   (b)  Formula 29 MATLAB code and its output:

    <>x0=2.3;
    <>h=0.1;
    <>d2=(exp(x0+2*h)-exp(x0-2*h))/(4*h)	 D2   =   10.0408
    <>d1=(exp(x0+h)-exp(x0-h))/(2*h)	 D1   =    9.9908
    <>fe=(4*d1-d2)/3, err=fe-exp(x0)
	 FE   =    9.974149167936696  ERR = -3.3287E-05
    <>h=0.01;
    <>d2=(exp(x0+2*h)-exp(x0-2*h))/(4*h)	 D2   =    9.9748
    <>d1=(exp(x0+h)-exp(x0-h))/(2*h)	 D1   =    9.9743
    <>fe=(4*d1-d2)/3, err=fe-exp(x0)
	 FE   =    9.974182451489680  ERR = -3.3250E-09

   (c) Using |f'''''(c)| < 12.2, we get |E_d| < h^4 * 12.2 / 30, which is
	evaluated by MATLAB as follows:
    <>h=0.1; h4=h*h*h*h; f5=12.2; ed=h4*f5/30	 ED    =    4.0667E-05
    <>h=0.01; h4=h*h*h*h; f5=12.2; ed=h4*f5/30	 ED    =    4.0667E-09


Exercise 16 of Section 6.1, p.338.

   (a) Formula 22 in MATLAB, and its output:

 <> num = [1 -8 0 8 -1]; dh=[-2 -1 0 1 2]; x0=3;

    h=0.1; x=x0*ones(num)+h*dh;
    f=[1.02962   1.06471   1.09861   1.13140   1.16315];
    long; df=(f*num')/(12*h), short; err= df-1/3
		 DF    =    .333324999999999	 ERR   =   -8.3333E-06

    h=0.001; x=x0*ones(num)+h*dh;
    f=[1.09795   1.09828   1.09861   1.09895   1.09928];
    long; df=(f*num')/(12*h), short; err= df-1/3
		 DF    =    .335833333333443	 ERR   =     .0025

   (b) Formula 24 can be estimated as follows:

	   E = E_r + E_d, so |E| <= |E_r| + |E_d|.

   Now for f(x)=ln(x), we have f''''(x) = 24/x^5, which is no bigger than
   24/(x0-2*h)^5 on the interval x-2h<c<x+2h, so we use that value to estimate
   |f''''(c)| when computing

	   |E_d|	<= 24*h^4/(30*(x0-2*h)^5).

   We are given that |e_k|<= 5*10^{-6}, for all k.  E_r from Formula 24 needs
   e_{-2}, e_{-1}, e_{1}, and e_{2}.  These can be estimated with the triangle
   inequality:

	   |E_r|	<= |-e_{-2}+8e_{-1}-8 e_{1}+e_{2}| /  (12*h)
		   <= |e_{-2}|+8|e_{-1}|+8|e_{1}|+|e_{2}| /  (12*h)
		   <= 18*5*10^{-6} /  (12*h)

   Thus |E| <= 24*h^4/(30*(x0-2*h)^5) + 18*5*10^{-6} / (12*h).  We compute
   this with MATLAB:

    <>h=0.1; h4=h*h*h*h; x0=3; x=x0-2*h; x5=x*x*x*x*x;
    <>eps=5e-6;   ed = 24*h4/(30*x5), er=18*eps/(12*h),  e=ed+er
    ED    =    4.648360802046769E-07
    ER    =    7.499999999999990E-05
    E     =    7.546483608020458E-05
    <>h=0.001; h4=h*h*h*h; x0=3; x=x0-2*h; x5=x*x*x*x*x;
    <>eps=5e-6;   ed = 24*h4/(30*x5), er=18*eps/(12*h),  e=ed+er
    ED    =    3.303176988919223E-15
    ER    =    .007500000000000
    E     =    .007500000000003



Algorithm 2(a,b,c,d,e) of Section 6.1, p.338.

  NOTE: there is a typo, relative error R[k]=2*E[k]/(|D[k]|+|D[k-1]|+eps)
  rather than 2*E[k]*(|D[k]|+|D[k-1]|+eps) in Program 6.1.

  We modify Program 6.1, p.332, to use the O(h^4) centered difference
  formula for the derivative.  In all cases, the iteration terminates
  because the total error starts to decrease at some not-so-small h value:

   (a)
   Total error has begun to increase.
   Iteration 6: h=1.0e-06, f'(x0)= 7.517349470731459e+01, E=1.3e-08, R=1.7e-10
   (b)
   Total error has begun to increase.
   Iteration 5: h=1.0e-05, f'(x0)= 1.228597423403680e+00, E=2.5e-11, R=2.0e-11
   (c)
   Total error has begun to increase.
   Iteration 6: h=1.0e-06, f'(x0)= 1.951559609023508e+00, E=1.0e-10, R=5.2e-11
   (d)
   Total error has begun to increase.
   Iteration 6: h=1.0e-06, f'(x0)= 2.965514829348740e+00, E=1.8e-10, R=6.0e-11
   (e)
   Total error has begun to increase.
   Iteration 9: h=1.0e-09, f'(x0)= 1.015206105898848e+00, E=1.8e-12, R=1.8e-12

  C PROGRAM:
   #include <stdio.h>
   #include <math.h>
   #define MAX1 15
   /* Set which part of the problem we will solve: */
   #define	X0	X0a;
   #define	f	fa;

   #define	X0e	0.0001;
   double	fe(double x){  return pow(x,pow(x,x));}

   #define	X0d	(1-sqrt(5.))/2;
   double	fd(double x){  return sin(pow(x,3.)-7*pow(x,2.)+6*x+8);}

   #define	X0c	1/sqrt(2.);
   double	fc(double x){  return sin(cos( 1/x));}

   #define	X0b	(1+sqrt(5.))/3;
   double	fb(double x){ return tan(cos((sqrt(5.)+sin(x))/(1+x*x)));}

   #define	X0a	(1/sqrt(3.))
   double	fa(double x){
     return 60*pow(x,45.)-32*pow(x,33.)+233*pow(x,5.)-47*x*x-77.;
   }

   void	message(const char *text)  { printf( "%s.\n", text); }

   /* O(h^4) derivative formula: */
   double   df(double x, double h)   {
     return ( -f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h) ) / ( 12*h );
   }

   int
   main(void)
   {
     double h, x, tolerance, eps, D[MAX1], H[MAX1], R[MAX1], E[MAX1];
     int k;

     eps = 1.0e-15;		/* Add this to avoid relative error problems */
     tolerance = 1.0e-13;	/* Target accuracy for the derivative */

     x=X0;			/* Point at which to differentiate */

     h=1;			/* Initial h */
     D[0]=df(x,h);		/* First approximation to the derivative */
     E[0]=R[0]=0;		/* Initialize error trackers */
     for(k=1; 1 ; k++) {
       h /= 10;			/* Decrease h for the next trial */
       D[k]=df(x,h);    E[k]=fabs(D[k]-D[k-1]);
       R[k]=2*E[k] / ( fabs(D[k]) + fabs(D[k-1]) + eps );
       if ( R[k]<tolerance ) {
	 message("Derivative found within desired tolerance"); break;    }
       if( k>1 && E[k]>E[k-1]) {
	 message("Total error has begun to increase"); break;    }
       if( k==MAX1-1  ) {
	 message("Maximum number of iterations reached"); break;    }
     }
     printf("Iteration %d: h=%6.1e, f'(x0)= %20.15e, E=%6.1e, R=%6.1e\n",
	    k, h, D[k], E[k], R[k] );
     return 0;
   }



Exercise 4 of Section 6.2, p.349.

   NOTE: part b has a typo: Should be "use h=0.1"  rather than "use h=0.01".

   MATLAB Code:

   (a)
    num=[1 -2 1]; x0=1; 
    h=0.05; short;  f=[.5817     .5403     .4976]
    long; d2f=(f*num')/(h*h), -cos(x0)
    short; err=d2f+cos(x0)
    D2F   =   -.520000000000009
    ANS   =   -.540302305868140		 ERR   =      .0203

   (b)
    num=[1 -2 1]; x0=1; 
    h=0.1; f=[.6216     .5403     .4536];
    long; d2f=(f*num')/(h*h), -cos(x0)
    short; err=d2f+cos(x0)
    D2F   =   -.539999999999996
    ANS   =   -.540302305868140		 ERR   =    3.0231E-04

   (c)
    num=[-1 16 -30 16 -1]; x0=1; 
    h=0.05; f=[.6216     .5817     .5403     .4976     .4536];
    long; d2f=(f*num')/(12*h*h), -cos(x0)
    short; err=d2f+cos(x0)
    D2F   =   -.513333333333358
    ANS   =   -.540302305868140		 ERR   =     .0270

   (c')
    num=[-1 16 -30 16 -1]; x0=1; 
    h=0.1; f=[.6967     .6216     .5403     .4536     .3624];
    long; d2f=(f*num')/(12*h*h), -cos(x0)
    short; err=d2f+cos(x0)
    D2F   =   -.540833333333323
    ANS   =   -.540302305868140		 ERR   =   -5.3103E-04

   (d) Value b is the most accurate of a,b,c.  Value c' is computed to show
     how the error actually decreases, even for the higher-order formula,
     when the step size h is increased.


Exercise 7 of Section 6.2, p.349.

  Subtracting:
  f(x+h)=f(x) + hf'(x) + h^2f''(x)/2 + h^3f'''(x)/6 + h^4f''''(x)/24 +O(h^5)
  f(x-h)=f(x) - hf'(x) + h^2f''(x)/2 - h^3f'''(x)/6 + h^4f''''(x)/24 +O(h^5)
  ---------------------------------------------------------------------------
  f(x+h)-f(x-h) = 2hf'(x) + h^3f'''(x)/3 +O(h^5)

  Subtracting:
  f(x+2h)=f(x) + 2hf'(x) + 2h^2f''(x) + 4h^3f'''(x)/3 + 2h^4f''''(x)/3 +O(h^5)
  f(x-2h)=f(x) - 2hf'(x) + 2h^2f''(x) - 4h^3f'''(x)/3 + 2h^4f''''(x)/3 +O(h^5)
  ---------------------------------------------------------------------------
  f(x+2h)-f(x-2h) = 4hf'(x) + 8h^3f'''(x)/3 +O(h^5)

  Eliminate f' by another subtraction:
  f(x+2h)-f(x-2h)  = 4hf'(x) + 8h^3f'''(x)/3 +O(h^5)
  2[f(x+h)-f(x-h)] = 4hf'(x) + 2h^3f'''(x)/3 +O(h^5)
  ---------------------------------------------------------------------------
  f(x+2h)-f(x-2h) - 2[f(x+h)-f(x-h)] = 6h^3f'''(x)/3 +O(h^5)

  Thus,
	 f'''(x) = [f(x+2h) - 2f(x+h) + 2f(x-h) - f(x-2h)]/2h^3 + O(h^2)



Algorithm 1 of Section 6.2, p.351.

  C Code:

  #include <assert.h>
  #include <math.h>
  #include <stdio.h>

  /* Differentiation formulas based on Newton polynomial interpolation. */

  /* Compute the Newton polynomial leading coefficients a[0],...,a[n] using
     the divided differences derived from {(x[k],f[k]) : k=0,1,...,n}: */
  void
  divided_differences ( double *a, const double *x, const double *f, int n)
  {
    int k, j;
    for(k=0; k<=n; k++) a[k] = f[k];  /* Initialize a[] with f's: */
    /* Compute jth divided difference, putting it in a[j]: */
    for(j=1; j<=n; j++)
      for(k=n; k>=j; k--)      a[k] = (a[k]-a[k-1])/(x[k]-x[k-j]);
    return;
  }
  /* Use Horner's method: */
  double
  evaluate_Newton_polynomial(const double *a, const double *x, int n, double t)
  {
    int k;  double y;
    y = a[n];
    for(k=n-1; k>=0; k--)    y = a[k] + (t-x[k])*y;
    return y;
  }
  /* Permute x0 <- x1 <- ... <- xn <- x0 */
  void
  cycle( double *z, int n)
  {
    double temp;  int k;
    temp = z[0];
    for(k=0; k<n; k++) z[k] = z[k+1];
    z[n] = temp;
    return;
  }
  /* Compute the Newton polynomial through collinear points and differentiate
     it at each of them: */
  #define N	7		/* Reasonably large */
  #define SIZE	N+1
  int
  main(void)
  {
    FILE *outfile;		/* Plot the Newton polynomial here */
    double a[SIZE]={0};		/* vector of Newton polynomial coefficients */
    double x[SIZE]={0, 1, 2, 3,  4,  5,  6,  7};
    double y[SIZE]={2, 4.2, 6,8,9.9,12,14.1,16}; /* collinear, with errors */
    int i;  double df, p, t, dt;

    /* Test if the Newton polynomial is correctly computed: */
    divided_differences(a,x,y,N); /* compute Newton polynomial */
    for(i=0; i<=N; i++) {
      p = evaluate_Newton_polynomial(a,x,N,x[i]);
      printf("P(x[%d])=P(%g)=%17.14f; y[%d]=%17.14f; |diff|=%7.1e\n",
	     i, x[i], p, i, y[i], fabs(p-y[i]) ); 
    }
    /* To get f'(x0), evaluate Q(x0), where Q is the Newton polynomial of
       degree n-1 determined by coefficients a[1],...,a[n] and centers
       x[1],...,x[n-1] */
    for(i=0; i<=N; i++) {
      divided_differences(a,x,y,N); /* compute Newton polynomial */
      df = evaluate_Newton_polynomial(a+1,x+1,N-1,x[0]); /* compute f'(x0) */
      printf(" f'(x[%d]) = f'(%g) = %17.14f\n", i, x[0], df );
      cycle(x,N);    cycle(y,N);  /* permute so that next x is in x[0]  */
    }
    /* Plot 100 values of p to illustrate polynomial wiggle: */
    outfile = fopen("6-2-a1.txt","w");  assert(outfile);
    printf("Plotting 100 points on the polynomial interpolating {(x,y)}.\n");
    dt = (x[N]-x[0])/100;
    for(i=0; i<100; i++) {
      t = x[0] + i*dt;
      fprintf(outfile, "%f %f\n",
	      t, evaluate_Newton_polynomial(a,x,N,t));
    }
    fclose(outfile);
    return 0;
  }

  OUTPUT:
     % cc -Wall 6-2-a1.c
     % ./a.out
     P(x[0])=P(0)= 2.00000000000000; y[0]= 2.00000000000000; |diff|=0.0e+00
     P(x[1])=P(1)= 4.20000000000000; y[1]= 4.20000000000000; |diff|=0.0e+00
     P(x[2])=P(2)= 6.00000000000000; y[2]= 6.00000000000000; |diff|=0.0e+00
     P(x[3])=P(3)= 8.00000000000000; y[3]= 8.00000000000000; |diff|=0.0e+00
     P(x[4])=P(4)= 9.90000000000000; y[4]= 9.90000000000000; |diff|=0.0e+00
     P(x[5])=P(5)=12.00000000000000; y[5]=12.00000000000000; |diff|=0.0e+00
     P(x[6])=P(6)=14.10000000000000; y[6]=14.10000000000000; |diff|=0.0e+00
     P(x[7])=P(7)=16.00000000000000; y[7]=16.00000000000000; |diff|=0.0e+00
      f'(x[0]) = f'(0) =  4.15833333333333
      f'(x[1]) = f'(1) =  1.56333333333333
      f'(x[2]) = f'(2) =  2.00833333333333
      f'(x[3]) = f'(3) =  1.92666666666667
      f'(x[4]) = f'(4) =  1.95166666666667
      f'(x[5]) = f'(5) =  2.21666666666667
      f'(x[6]) = f'(6) =  1.85500000000000
      f'(x[7]) = f'(7) =  2.70000000000000
     Plotting 100 points on the polynomial interpolating {(x,y)}.