roninkoi

Bisection and Newton roots

Feb 22nd, 2021
1,011
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <math.h>
  4.  
  5. // sign of floating point number
  6. #define sign(x) (x >= 0.0 ? 1.0 : -1.0)
  7.  
  8. // function f
  9. #define f(x) (sin(3 * M_PI * (x) * (x) * (x) / ((x) * (x) - 1.0)) + 0.5)
  10.  
  11. // Bisection method tolerance
  12. #define BISECT_TOL 1.0e-7
  13. // Bisection method maximum iterations
  14. #define BISECT_MAXIT 10000
  15.  
  16. /*
  17.  * Calculates root of equation f(x) using bisection method
  18.  * in the interval [a, b].
  19.  */
  20. double bisect_f(double a, double b)
  21. {
  22.     double c; // solution to f(c) == 0
  23.     double aa = a; // current interval
  24.     double bb = b;
  25.  
  26.     int conv = 0;
  27.    
  28.     for (int i = 0; i < BISECT_MAXIT; ++i) {
  29.         c = (aa + bb) / 2.0; // midpoint of current interval
  30.  
  31.         double fc = f(c);
  32.  
  33.         if (fc == 0.0 || (b - a) / 2.0 < BISECT_TOL) {
  34.             conv = 1; // found root!
  35.             break;
  36.         }
  37.  
  38.         if (sign(fc) == sign(f(aa))) // pick which side to advance
  39.             aa = c;
  40.         else
  41.             bb = c;
  42.     }
  43.  
  44.     if (!conv)
  45.         fprintf(stderr, "Bisect error: Didn't find root!\n");
  46.  
  47.     return c;
  48. }
  49.  
  50. // Newton's method maximum iterations
  51. #define NEWTON_MAXIT 10000
  52. // Newton's method tolerance
  53. #define NEWTON_TOL 1.0e-7
  54. // Newton's method smallest value
  55. #define NEWTON_EPS 1.0e-10
  56. // Newton's method derivative dx
  57. #define NEWTON_DX 1.0e-7
  58.  
  59. /*
  60.  * Calculates root of f(x) using Newton's method
  61.  * based on initial guess x0. Approximates derivative f'(x).
  62.  */
  63. double newton_f(double x0)
  64. {
  65.     double x = x0;
  66.  
  67.     int conv = 0;
  68.     for (int i = 0; i < NEWTON_MAXIT; ++i) {
  69.         double y = f(x);
  70.         double dy = (f(x + NEWTON_DX) - f(x - NEWTON_DX)) / NEWTON_DX / 2.0; // approximate derivative
  71.  
  72.         if (fabs(dy) <= NEWTON_EPS) {
  73.             fprintf(stderr, "Newton error: Small f'(x)!\n");
  74.             break;
  75.         }
  76.  
  77.         double x1 = x - y / dy; // advance
  78.  
  79.         if (fabs(x1 - x) <= NEWTON_TOL) {
  80.             conv = 1; // converged
  81.             break;
  82.         }
  83.  
  84.         x = x1; // new value of x
  85.     }
  86.  
  87.     if (!conv)
  88.         fprintf(stderr, "Newton error: Didn't converge!\n");
  89.  
  90.     return x;
  91. }
  92.  
RAW Paste Data