SWMM5 findroot.c, Summary

 Below is a step‐by‐step explanation of findroot.c, which provides two root-finding methods for solving 

$f(x) = 0$: a Newton‐Raphson method (with fallback bisection) and Ridder’s Method. Both functions iterate until either the root is found within a given accuracy or a maximum iteration limit is reached.

---

1. Overall Purpose

//  findroot.c
//
//  Finds solution of func(x) = 0 using either ...

• This file is for finding roots of functions $f(x)$ = 0, i.e., we want to find $x$ such that $f(x)$ is zero within a tolerance.

• It contains two main functions:

  1. findroot_Newton(...): A combination of Newton–Raphson and bisection.
  2. findroot_Ridder(...): Implementation of Ridder’s Method.

• The user supplies:

  • A bracket (x1, x2) such that $f(x1)$ and $f(x2)$ have opposite signs (i.e., the function’s sign changes, indicating a root in between).
  • A function pointer that either returns $f(x)$ & derivative or just $f(x)$.

---

2. findroot_Newton(...)

int findroot_Newton(
    double x1, double x2, double* rts, double xacc,
    void (*func) (double x, double* f, double* df, void* p),
    void* p)

1. Parameters:

   • x1, x2: The initial bracket endpoints.
   • *rts: On input, an initial guess for the root; on output, the refined root.
   • xacc: The required accuracy for the solution.
   • func(x, &f, &df, p): A user function that returns:
     • f: the function value $f(x)$,
     • df: the derivative $f'(x)$.
   • p: A pointer to any extra data needed by func(). (It can be NULL if no extra data is needed.)

2. Algorithm:

   1. Initial Setup:
      • xlo = x1, xhi = x2, x = *rts.
      • Evaluate func(x, &f, &df, p).
      • Keep track of iteration count in n.
   2. Loop (up to MAXIT=60 times):
      • The method tries a Newton‐Raphson step: dx = f/df.
      • But checks if that step is out of range or not decreasing well:

        if ( ( (x - xhi)*df - f)*((x - xlo)*df - f) >= 0.0
          || (fabs(2.0*f) > fabs(dxold*df)) )
        {
          // If so, do a bisection step
          dx = 0.5*(xhi - xlo);
          x = xlo + dx;
        }
        else
        {
          // Otherwise, do Newton step
          dx = f/df;
          x -= dx;
        }
        

      • After adjusting x, check if |dx| < xacc (i.e., the change is within tolerance).
      • Evaluate func(x, &f, &df, p) again.
      • If f < 0.0, set xlo = x; else xhi = x.
   3. Convergence:
      • If it converges or if maximum iterations are reached, store x in *rts.
      • The function returns the number of function evaluations used or 0 if MAXIT was exceeded (i.e., no convergence in time).

3. Notes:

   • This approach tries Newton’s method if possible (fast convergence) but falls back to bisection if the Newton step might overshoot or not reduce the bracket effectively.

---

3. findroot_Ridder(...)

double findroot_Ridder(
    double x1, double x2, double xacc,
    double (*func)(double, void* p), void* p)

1. Parameters:

   • x1, x2: The bracket endpoints (sign of f(x1) and f(x2) must differ).
   • xacc: Accuracy required.
   • func(x, p): The function returning only f(x).
   • p: Optional user data pointer.

2. Algorithm: Ridder’s method is a bracketing method that uses intermediate function evaluations and a special formula to find a better midpoint.

   1. Evaluate flo = func(x1, p), fhi = func(x2, p).
   2. If flo == 0.0, root is x1; if fhi == 0.0, root is x2.
   3. Check (flo > 0.0 && fhi < 0.0) || (flo < 0.0 && fhi > 0.0) to confirm it’s bracketed.
   4. Then iterate up to MAXIT=60:
      • xm = midpoint: (xlo + xhi)/2.
      • Evaluate fm = func(xm, p).
      • Compute s = sqrt( fm^2 - flo * fhi ).
      • If s == 0.0 => can’t proceed, return the midpoint.
      • Then a new approximation xnew is found by a Ridder formula: $$xnew = xm + (xm - xlo) * \bigl( (\text{flo} \ge \text{fhi} ? 1 : -1)*fm / s \bigr)$$
      • If |xnew - ans| <= xacc, we break (converged).
      • Evaluate fnew = func(xnew,p).
      • Then we decide how to re-bracket:

        if ( SIGN(fm, fnew) != fm ) {
            xlo = xm; flo = fm;
            xhi = xnew; fhi = fnew;
        }
        else if ( SIGN(flo, fnew) != flo ) {
            xhi = xnew; fhi = fnew;
        }
        else if ( SIGN(fhi, fnew) != fhi ) {
            xlo = xnew; flo = fnew;
        }
        

      • If bracket becomes very small or we converge, we return.
   5. If no bracket or no convergence, returns -1.e20 indicating failure.

3. Notes:

   • Ridder’s method is a robust bracketing approach with guaranteed convergence if sign change is present, typically faster than simple bisection but less complicated than Newton.

---

4. Implementation Details

• MAXIT = 60 is the maximum iteration count for both methods.
• SIGN(a,b) macro returns |a| if b >= 0, else -|a|.
• Each function uses the bracket [x1, x2] and ensures it always encloses the root (i.e., f(xlo)*f(xhi) < 0).
• findroot_Newton(...) uses a “mixed” approach: Newton’s method but does a bisection if the Newton step is questionable.
• findroot_Ridder(...) purely uses Ridder’s formula, a special bracket-based improvement over bisection.

---

5. Key Takeaways

• These functions assume the user has already found an initial bracket [x1, x2] with f(x1)*f(x2) < 0.
• findroot_Newton() returns an integer: the count of function evaluations or zero if it fails.
• findroot_Ridder() returns a double root approximation or -1.e20 if bracket was invalid or no success.
• Both methods can be integrated easily by providing a function pointer for f(x) (and derivative if needed).

Hence, findroot.c is a utility file for any part of the SWMM code (or related code) needing a robust root-finding procedure.