Solutions to Homework Assignment 2

Here are my solutions to homework 2. Please let me know if you find an error!

Problem 7.8 in Austin/Chancogne.

/*
 *  ======================================================================
 *  prog_cuberoot.c : Compute cube root of a number           
 *
 *  Written By : Rachel Albrecht and Mark Austin              October 1996
 *  ======================================================================
 */

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

int main(void){
float fN, fCr;
float CubeRoot(float fN);

   /* [a] : prompt user for number */

   printf("\n");
   printf("Please enter a value of N : ");
   scanf("%f%*c",&fN);
   printf("The value of N is : %f\n",fN);

   /* [b] : compute and print cube root */

   fCr = CubeRoot (fN);
   printf("The cube root of N is : %f\n",fCr );
}

/*
 *  ======================================================================
 *  CubeRoot() : Compute cube root of a number           
 *
 *  Input  : float fPassed -- Number that cube root is to be computed.
 *  Output : float fX      -- Cube root of fPassed.
 *  ======================================================================
 */

#define EPSILON 0.00001

float CubeRoot(float fPassed) {
int iCount = 0, iConvergenceAchieved;
float fX, fXnext;

   /* [a] : No iteration is required for trivial cases */

   if (fPassed == 1 || fPassed == 0 || fPassed == -1 ){
       return fPassed;
   }

   /* [b] : Compute cube root */

   fX = fPassed/2.0;  /* initial guess for cube root of fPassed */

   iConvergenceAchieved = 0;
   while(iConvergenceAchieved == 0) {

        /* [b.1] : Compute next estimate of cube root */

        fXnext = (2*fX*fX*fX + fPassed)/(3*fX*fX);

        /* [b.3] : Has algorithm converged */

        if(fabs((fXnext -fX)/fX) < EPSILON )
           iConvergenceAchieved = 1;

        /* [b.3] : Update variables */

        fX     = fXnext;
        iCount = iCount + 1;
   }

   /* [c] : Print results on convergence and return result */

   printf("Cube Root algorithm converged in %3d iterations\n", iCount );

   return fX;
}

Example. Here is a typical script of program I/O;

     prompt >>
     prompt >> CUBEROOT
     Please enter a value of N : 20
     The value of N is : 20.000000
     Cube Root algorithm converged in 7 iterations
     The cube root of N is : 2.714417
     prompt >>

Problem 6.3 in Austin/Chancogne.

/*
 *  ========================================================================
 *  prog_sphere.c -- Use method of Newton-Raphson to compute the depth at
 *                   a sphere of density "rho" sits in water.
 *
 *  Written By : Mark Austin                                    October 1996 
 *  ========================================================================
 */

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

double NewtonRaphson( double, double, double, double, int );
double dDensity;  /* Global variable for density of sphere */

int main( void ) {
enum { NoPoints = 11 };
double dResponse[NoPoints][2];
double dRoot;
float fStart;
int ii;

   /* [a] : Compute depth of sphere for rho = 0 to rho = 1 */

   for(ii = 1; ii <= NoPoints; ii = ii + 1) {

       /* [a.1] : Empirically adjust starting value for iteration */

       if(ii <= 9)
         fStart = 1.0;
       else 
         fStart = 2.5;

       /* [a.2] : Compute root to equation of equilibrium */

       dDensity = (ii-1)*0.1;
       dRoot = NewtonRaphson( fStart, 0.00001, 0.00001, 0.1, 20 );

       /* [a.3] : Save density and computed depth */

       dResponse[ii-1][0] = dDensity;
       dResponse[ii-1][1] = dRoot; 
   }

   /* [b] : Print results of analysis */

   printf("==================\n");
   printf("Density      [d/r]\n");
   printf("==================\n");

   for(ii = 1; ii <= NoPoints; ii = ii + 1) 
       printf("%7.2f  %9.2f\n", dResponse[ii-1][0], dResponse[ii-1][1]);
} 

/*
 *  =====================================================
 *  Function for equilibrium of sphere floating in water.
 *  =====================================================
 */

double SphereEquilibrium(double dX) {
       return (dX*dX*dX  - 3*dX*dX + 4.0*dDensity);
}

/*
 *  ========================================================================
 *  Use Newton Raphson to Compute Root of Equation
 *
 *  Note : First and second derivatives are estimated via finite differences
 *       : We use the second derivative in test for convergence criteria
 *
 *  Input  : dX0        -- Initial Guess at Root.
 *         : dDelta     -- Convergence criteria for change in x values
 *         : dEpsilon   -- Convergence criteria for y-values.
 *         : dIncrement -- Increment for derivative estimate
 *         : iMaxIter   -- Maximum number of iteration
 *  Output : dXnew      -- Root of the equation   
 *  ========================================================================
 */

double NewtonRaphson( double dX0, double dDelta, double dEpsilon,
                      double dIncrement, int iMaxIter) {

double dXn, dXnew, dYn, dAbsXn, dAbsXndiff, dAbsYn;
double d1Derivative, d2Derivative;
int iCount = 0;

   /* [a] : Initialize variables for iteration */

   dXn        = dX0;
   dAbsXn     = fabs(dXn);
   dAbsXndiff = fabs(dXn - dXnew);
   dAbsYn     = fabs( SphereEquilibrium(dXn) );

   /* [b] : Compute Newton-Raphson Iterates */
 
   while((dAbsYn > dEpsilon) && (iCount <= iMaxIter) &&
         (dAbsXndiff > dDelta*dAbsXn) ) {

      /* [b.1] : Compute finite difference estimates of first and  */
      /*         second derivates                                  */

      d1Derivative = ( SphereEquilibrium(dXn + dIncrement) - 
                       SphereEquilibrium(dXn - dIncrement) )/(2.0*dIncrement);
      d2Derivative = ( SphereEquilibrium(dXn + dIncrement) -
                       2.0*SphereEquilibrium (dXn)
                     + SphereEquilibrium(dXn - dIncrement) )/pow(dIncrement, 2.0);

      /* [b.2] : If the sufficient condition for the algorithm to convergence */
      /*         is not satisfied, send error message and exit program.       */

      if ( fabs( SphereEquilibrium(dXn))*d2Derivative >= pow(d1Derivative, 2.0)) {
         printf("ERROR >> Algorithm may not converge for the starting value.\n");
         printf("ERROR >> Please choose another starting value.\n\n");
         exit(1);
      }  

      /* [b.3] : Compute update in Newton-Raphson iteration */

      dXnew  = dXn - SphereEquilibrium(dXn)/d1Derivative;
      dAbsXn = fabs(dXn);
      dAbsXndiff = fabs(dXnew - dXn);
      dAbsYn     = fabs( SphereEquilibrium (dXn) );
        
      iCount += 1;
      dXn = dXnew; 
   }

   return dXn;
}

The function SphereEquilibrium() computes the equilibrium of the sphere floating in water. Calls to SphereEquilibrium() are made from {\tt NewtonRaphson()} as needed. A more elegant solution -- beyond the scope of this book, unfortunately -- would use pointers to functions within ewtonRaphson() and pass the address of SphereEquilibrium()} as an argument to NewtonRaphson().

Example. This program generates the output:

    ==================
    Density      [d/r]
    ==================
       0.00      -0.00
       0.10       0.39
       0.20       0.57
       0.30       0.73
       0.40       0.87
       0.50       1.00
       0.60       1.13
       0.70       1.27
       0.80       1.43
       0.90       2.35
       1.00       2.00