findroot.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // findroot.c | ||
| 3 | // | ||
| 4 | // Finds solution of func(x) = 0 using either the Newton-Raphson | ||
| 5 | // method or Ridder's Method. | ||
| 6 | // Based on code from Numerical Recipes in C (Cambridge University | ||
| 7 | // Press, 1992). | ||
| 8 | // | ||
| 9 | // Date: 11/19/13 | ||
| 10 | // Author: L. Rossman | ||
| 11 | //----------------------------------------------------------------------------- | ||
| 12 | |||
| 13 | #include <math.h> | ||
| 14 | #include "findroot.h" | ||
| 15 | |||
| 16 | #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) | ||
| 17 | #define MAXIT 60 | ||
| 18 | |||
| 19 | |||
| 20 | 1429422 | int findroot_Newton(double x1, double x2, double* rts, double xacc, | |
| 21 | void (*func) (double x, double* f, double* df, void* p), | ||
| 22 | void* p) | ||
| 23 | // | ||
| 24 | // Using a combination of Newton-Raphson and bisection, find the root of a | ||
| 25 | // function func bracketed between x1 and x2. The root, returned in rts, | ||
| 26 | // will be refined until its accuracy is known within +/-xacc. func is a | ||
| 27 | // user-supplied routine, that returns both the function value and the first | ||
| 28 | // derivative of the function. p is a pointer to any auxilary data structure | ||
| 29 | // that func may require. It can be NULL if not needed. The function returns | ||
| 30 | // the number of function evaluations used or 0 if the maximum allowed | ||
| 31 | // iterations were exceeded. | ||
| 32 | // | ||
| 33 | // NOTES: | ||
| 34 | // 1. The calling program must insure that the signs of func(x1) and func(x2) | ||
| 35 | // are not the same, otherwise x1 and x2 do not bracket the root. | ||
| 36 | // 2. If func(x1) > func(x2) then the order of x1 and x2 should be | ||
| 37 | // switched in the call to Newton. | ||
| 38 | // | ||
| 39 | { | ||
| 40 | 1429422 | int j, n = 0; | |
| 41 | double df, dx, dxold, f, x; | ||
| 42 | double temp, xhi, xlo; | ||
| 43 | |||
| 44 | // Initialize the "stepsize before last" and the last step. | ||
| 45 | 1429422 | x = *rts; | |
| 46 | 1429422 | xlo = x1; | |
| 47 | 1429422 | xhi = x2; | |
| 48 | 1429422 | dxold = fabs(x2-x1); | |
| 49 | 1429422 | dx = dxold; | |
| 50 | 1429422 | func(x, &f, &df, p); | |
| 51 | 1429422 | n++; | |
| 52 | |||
| 53 | // Loop over allowed iterations. | ||
| 54 |
1/2✓ Branch 0 taken 7243254 times.
✗ Branch 1 not taken.
|
7243254 | for (j=1; j<=MAXIT; j++) |
| 55 | { | ||
| 56 | // Bisect if Newton out of range or not decreasing fast enough. | ||
| 57 |
2/2✓ Branch 0 taken 6311799 times.
✓ Branch 1 taken 931455 times.
|
7243254 | if ( ( ( (x-xhi)*df-f)*((x-xlo)*df-f) >= 0.0 |
| 58 |
2/2✓ Branch 0 taken 1020 times.
✓ Branch 1 taken 6310779 times.
|
6311799 | || (fabs(2.0*f) > fabs(dxold*df) ) ) ) |
| 59 | { | ||
| 60 | 932475 | dxold = dx; | |
| 61 | 932475 | dx = 0.5*(xhi-xlo); | |
| 62 | 932475 | x = xlo + dx; | |
| 63 |
2/2✓ Branch 0 taken 249 times.
✓ Branch 1 taken 932226 times.
|
932475 | if ( xlo == x ) break; |
| 64 | } | ||
| 65 | |||
| 66 | // Newton step acceptable. Take it. | ||
| 67 | else | ||
| 68 | { | ||
| 69 | 6310779 | dxold = dx; | |
| 70 | 6310779 | dx = f/df; | |
| 71 | 6310779 | temp = x; | |
| 72 | 6310779 | x -= dx; | |
| 73 |
2/2✓ Branch 0 taken 13755 times.
✓ Branch 1 taken 6297024 times.
|
6310779 | if ( temp == x ) break; |
| 74 | } | ||
| 75 | |||
| 76 | // Convergence criterion. | ||
| 77 |
2/2✓ Branch 0 taken 1415418 times.
✓ Branch 1 taken 5813832 times.
|
7229250 | if ( fabs(dx) < xacc ) break; |
| 78 | |||
| 79 | // Evaluate function. Maintain bracket on the root. | ||
| 80 | 5813832 | func(x, &f, &df, p); | |
| 81 | 5813832 | n++; | |
| 82 |
2/2✓ Branch 0 taken 871769 times.
✓ Branch 1 taken 4942063 times.
|
5813832 | if ( f < 0.0 ) xlo = x; |
| 83 | 4942063 | else xhi = x; | |
| 84 | } | ||
| 85 | 1429422 | *rts = x; | |
| 86 |
1/2✓ Branch 0 taken 1429422 times.
✗ Branch 1 not taken.
|
1429422 | if ( n <= MAXIT) return n; |
| 87 | ✗ | else return 0; | |
| 88 | }; | ||
| 89 | |||
| 90 | |||
| 91 | 329238 | double findroot_Ridder(double x1, double x2, double xacc, | |
| 92 | double (*func)(double, void* p), void* p) | ||
| 93 | { | ||
| 94 | int j; | ||
| 95 | double ans, fhi, flo, fm, fnew, s, xhi, xlo, xm, xnew; | ||
| 96 | |||
| 97 | 329238 | flo = func(x1, p); | |
| 98 | 329238 | fhi = func(x2, p); | |
| 99 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 329238 times.
|
329238 | if ( flo == 0.0 ) return x1; |
| 100 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 329238 times.
|
329238 | if ( fhi == 0.0 ) return x2; |
| 101 | 329238 | ans = 0.5*(x1+x2); | |
| 102 |
7/8✓ Branch 0 taken 27617 times.
✓ Branch 1 taken 301621 times.
✓ Branch 2 taken 288 times.
✓ Branch 3 taken 27329 times.
✓ Branch 4 taken 301621 times.
✓ Branch 5 taken 288 times.
✓ Branch 6 taken 301621 times.
✗ Branch 7 not taken.
|
329238 | if ( (flo > 0.0 && fhi < 0.0) || (flo < 0.0 && fhi > 0.0) ) |
| 103 | { | ||
| 104 | 328950 | xlo = x1; | |
| 105 | 328950 | xhi = x2; | |
| 106 |
1/2✓ Branch 0 taken 936038 times.
✗ Branch 1 not taken.
|
936038 | for (j=1; j<=MAXIT; j++) { |
| 107 | 936038 | xm = 0.5*(xlo + xhi); | |
| 108 | 936038 | fm = func(xm, p); | |
| 109 | 936038 | s = sqrt( fm*fm - flo*fhi ); | |
| 110 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 936038 times.
|
936038 | if (s == 0.0) return ans; |
| 111 |
2/2✓ Branch 0 taken 577110 times.
✓ Branch 1 taken 358928 times.
|
936038 | xnew = xm + (xm-xlo)*( (flo >= fhi ? 1.0 : -1.0)*fm/s ); |
| 112 |
2/2✓ Branch 0 taken 328950 times.
✓ Branch 1 taken 607088 times.
|
936038 | if ( fabs(xnew - ans) <= xacc ) break; |
| 113 | 607088 | ans = xnew; | |
| 114 | 607088 | fnew = func(ans, p); | |
| 115 |
4/4✓ Branch 0 taken 85617 times.
✓ Branch 1 taken 521471 times.
✓ Branch 2 taken 561636 times.
✓ Branch 3 taken 45452 times.
|
607088 | if ( SIGN(fm, fnew) != fm) |
| 116 | { | ||
| 117 | 561636 | xlo = xm; | |
| 118 | 561636 | flo = fm; | |
| 119 | 561636 | xhi = ans; | |
| 120 | 561636 | fhi = fnew; | |
| 121 | } | ||
| 122 |
4/4✓ Branch 0 taken 28310 times.
✓ Branch 1 taken 17142 times.
✓ Branch 2 taken 17142 times.
✓ Branch 3 taken 28310 times.
|
45452 | else if ( SIGN(flo, fnew) != flo ) |
| 123 | { | ||
| 124 | 17142 | xhi = ans; | |
| 125 | 17142 | fhi = fnew; | |
| 126 | } | ||
| 127 |
2/4✓ Branch 0 taken 28310 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 28310 times.
✗ Branch 3 not taken.
|
28310 | else if ( SIGN(fhi, fnew) != fhi) |
| 128 | { | ||
| 129 | 28310 | xlo = ans; | |
| 130 | 28310 | flo = fnew; | |
| 131 | } | ||
| 132 | ✗ | else return ans; | |
| 133 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 607088 times.
|
607088 | if ( fabs(xhi - xlo) <= xacc ) return ans; |
| 134 | } | ||
| 135 | 328950 | return ans; | |
| 136 | } | ||
| 137 | 288 | return -1.e20; | |
| 138 | } | ||
| 139 |