Functions as parameters

1. Functions as parameters Figure 12.4 (assume this is stored in file methods.h) #include <stdio.h> #include <math.h> #define TRACE #define FALSE 0 #define TRUE 1 #define MAX_APPROX 100 /*iterations allowed before divergence assumed */ /* * Implements Newton's method for finding a root of a function f. * Finds a root (and sets output parameter error flag to FALSE) if two * successive approximations are found that differ by less than epsilon. * Sets output parameter error flag to TRUE if fderiv (the derivative * of f) returns a value of zero or if method fails to converge in * MAX_APPROX iterations */ double newton(double f(double farg), /* input - the function */ double fderiv(double fdarg), /* input - the function's first derivative */ double guess, /* input - user's guess for a root */ double epsilon, /* input - error tolerance */ int *errp) /* output - error flag */ TDBA66 Lecture Ch 12 vt-03

2. Fig. 12.4 cont. { double previous, /* previous approximation */ approx, /* new approximation */ deriv_val; /* value of derivative */ int num_approx = 0; /* number of iterations so far */ *errp = FALSE; approx = guess; do { ++num_approx; previous = approx; deriv_val = fderiv(previous); if (num_approx > MAX_APPROX) *errp = TRUE; else if (deriv_val == 0) *errp = TRUE; else approx = previous - f(previous) / deriv_val; /* Numbers and displays new approximation if trace is on */ #if defined(TRACE) if (*errp) printf("Approximation does not converge\n"); else printf("%3d. Function value at %.7f = %e\n", num_approx, approx, f(approx)); #endif } while (!(*errp) && fabs(approx - previous) >= epsilon); return (approx); } TDBA66 Lecture Ch 12 vt-03

3. Fig. 12.3 How to call /* * Finds roots of the equations * f(x) = 0 and g(x) = 0 * given a user-supplied guess and applying Newton's method * from a personal numerical methods library. */ #include <stdio.h> #include <math.h> #include "methods.h" /* Header file for personal numerical methods library containing function newton */ double f(double x); double fprime(double x); double g(double x); double gprime(double x); double newton(double f(double farg, double fderiv(double fdarg), double guess, double epsilon, int *errp); int main(void) { double guess, /* user's guess for function root */ epsilon, /* error tolerance */ root; int error; /* Get guess and error tolerance from user for f root approximation */ printf("\nEnter guess for root of f> "); scanf("%lf", &guess); printf("\nEnter tolerance> "); scanf("%lf", &epsilon); TDBA66 Lecture Ch 12 vt-03

4. Fig. 12.3 cont. /* Use function newton to look for root of f */ printf("\nFunction f\n"); root = newton(f, fprime, guess, epsilon, &error); if (error) printf("Root of f not found for guess of %.7f\n", guess); else printf(" f(%.7f) = %e\n", root, f(root)); return 0; } /* Functions for which roots are sought and functions' derivatives */ /* 3 2 * 5x - 2x + 3 */ double f(double x) { return (5 * pow(x, 3.0) - 2 * pow(x, 2.0) + 3); } /* 2 * 15x - 4x */ double fprime(double x) { return (15 * pow(x, 2.0) - 4 * x); } TDBA66 Lecture Ch 12 vt-03

5. Figure 12.19 Euler Method Program (extended) /* * Euler method approximation of solution to initial value problem * * dy * ---- = -cy , (x , y ) = (0, 1) * dx * * with step size H, output interval = 1 (every STEPS_PER_YEAR steps) * * Traces decay of a radioactive material, printing fraction of * original radioactivity remaining for each year. */ #include <stdio.h> #include <math.h> #define H 0.01 /* step size */ #define STEPS_PER_YEAR 100 /* integer representing 1/H */ int main(void) { double halflife; /* input - half-life of the element in question */ int years; /* input - number of years user */ FILE *outfile; /* requests */ double y, /* output - approximate result */ exact_y, /* output - actual function value */ error, /* output - difference between actual and approximate values */ c; /* ln(2) / half-life */ int x, /* current year */ i; /* loop counter */ TDBA66 Lecture Ch 12 vt-03

6. Fig. 12.19 cont. /* Get half-life from user and compute c */ printf("Half-life in years?> "); scanf("%lf", &halflife); c = log(2) / halflife; /* Get number of years to simulate */ printf("How many years should I simulate?> "); scanf("%d", &years); /* Display table heading and first line, which contains initial value */ printf("\nTIME(YEARS) ACTUAL APPROX "); printf("ERROR(ACTUAL-APPROX)\n"); x = 0; y = 1; /* Initial values */ printf("%5d%18.6f\n", x, y); /* Display approximations and exact values for number of years requested. Also print on outfile (plot data) */ outfile=fopen(“plot_data.txt”, “w”); for (x = 1; x <= years; ++x) { /* Compute one year's approximation */ for (i = 0; i < STEPS_PER_YEAR; ++i) y += H * -c * y; /* Compute exact value and error for comparison */ exact_y = exp(-c * x); error = exact_y - y; printf("%5d%18.6f%12.6f%19.6e\n", x, exact_y, y, error); fprintf(outfile,"%12.6f%12.6f\n", x, y); } return (0); } TDBA66 Lecture Ch 12 vt-03

7. Run Half-life in years?> 5.26 How many years should I simulate?> 12 TIME(YEARS) ACTUAL APPROX ERROR(ACTUAL-APPROX) 0 1.000000 1 0.876536 0.876460 7.616971e-05 2 0.768316 0.768183 1.335252e-04 3 0.673457 0.673281 1.755520e-04 4 0.590310 0.590104 2.051613e-04 5 0.517428 0.517203 2.247795e-04 6 0.453544 0.453308 2.364226e-04 7 0.397548 0.397306 2.417614e-04 8 0.348465 0.348223 2.421753e-04 9 0.305443 0.305204 2.387996e-04 10 0.267732 0.267499 2.325638e-04 11 0.234676 0.234452 2.242260e-04 12 0.205702 0.205488 2.144004e-04 How to plot? Use gnuplot TDBA66 Lecture Ch 12 vt-03