Below is a step‐by‐step explanation of findroot.c, which provides two root-finding methods for solving
: 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 = 0, i.e., we want to find such that is zero within a tolerance.
-
It contains two main functions:
findroot_Newton(...): A combination of Newton–Raphson and bisection.findroot_Ridder(...): Implementation of Ridder’s Method.
-
The user supplies:
- A bracket (
x1,x2) such that and have opposite signs (i.e., the function’s sign changes, indicating a root in between). - A function pointer that either returns & derivative or just .
- A bracket (
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)
-
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 ,df: the derivative .
p: A pointer to any extra data needed byfunc(). (It can be NULL if no extra data is needed.)
-
Algorithm:
- Initial Setup:
xlo = x1,xhi = x2,x = *rts.- Evaluate
func(x, &f, &df, p). - Keep track of iteration count in
n.
- Loop (up to
MAXIT=60times):- 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, setxlo = x; elsexhi = x.
- The method tries a Newton‐Raphson step:
- Convergence:
- If it converges or if maximum iterations are reached, store
xin*rts. - The function returns the number of function evaluations used or 0 if
MAXITwas exceeded (i.e., no convergence in time).
- If it converges or if maximum iterations are reached, store
- Initial Setup:
-
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)
-
Parameters:
x1, x2: The bracket endpoints (sign off(x1)andf(x2)must differ).xacc: Accuracy required.func(x, p): The function returning onlyf(x).p: Optional user data pointer.
-
Algorithm: Ridder’s method is a bracketing method that uses intermediate function evaluations and a special formula to find a better midpoint.
- Evaluate
flo = func(x1, p),fhi = func(x2, p). - If
flo == 0.0, root isx1; iffhi == 0.0, root isx2. - Check
(flo > 0.0 && fhi < 0.0) || (flo < 0.0 && fhi > 0.0)to confirm it’s bracketed. - 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
xnewis found by a Ridder formula: - 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.
- If no bracket or no convergence, returns
-1.e20indicating failure.
- Evaluate
-
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|ifb >= 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]withf(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.e20if 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.
