zeroNewton.c File Reference

Detailed Description

Implements the zeroNewton function for finding zeros of a function using Newton's method with numerical derivatives.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, C. Saunders, R. Soliday

Definition in file zeroNewton.c.

#include "mdb.h"

Go to the source code of this file.

Functions
double	zeroNewton (double(*fn)(), double value, double x_i, double dx, long n_passes, double _zero)
	Finds the zero of a function using Newton's method with numerical derivative computation.

Function Documentation

zeroNewton()

double zeroNewton	(	double(*	fn )(),
		double	value,
		double	x_i,
		double	dx,
		long	n_passes,
		double	_zero )

Finds the zero of a function using Newton's method with numerical derivative computation.

This function attempts to find a value x such that fn(x) is approximately equal to the given value. It employs Newton's iterative method, where the derivative is estimated numerically using finite differences.

Parameters

fn	Pointer to the function for which the zero is to be found.
value	The target value to solve for, i.e., find x such that fn(x) = value.
x_i	Initial guess for the independent variable.
dx	Increment used for numerical derivative computation.
n_passes	Number of iterations to perform.
_zero	Acceptable tolerance for the zero, determining when to stop the iterations.

Returns

The x value where fn(x) is approximately equal to value, or DBL_MAX if no zero is found within the iteration limit.

Definition at line 33 of file zeroNewton.c.

  {
  double f1, f2;
  double x1, x2;
  double dfdx;
  long i;
 
  x1 = x2 = x_i;
  for (i = 0; i < n_passes; i++) {
    if (fabs(f1 = (*fn)(x1)-value) < _zero)
      return (x1);
    if (i == n_passes - 1)
      return ((x1 + x2) / 2);
    f2 = (*fn)(x2 = x1 + dx) - value;
    dfdx = ((f2 - f1) / dx);
    x1 += -f1 / dfdx;
  }
  return (DBL_MAX); /* shouldn't happen */
}