zeroInterp.c File Reference

Implements the zeroInterp function for finding zeros of a function using successive interpolation. More...

#include <math.h>
#include <stdio.h>
#include "mdb.h"

Go to the source code of this file.

Macros
#define	sign(x)

Functions
double	zeroInterp (double(*fn)(), double value, double x_i, double x_f, double dx, double _zero)
	Finds the zero of a function within a specified interval using successive interpolation.

Detailed Description

Implements the zeroInterp function for finding zeros of a function using successive interpolation.

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 zeroInterp.c.

Macro Definition Documentation

sign

#define sign

(

x

)

Value:

((x) > 0 ? 1 : ((x) == 0 ? 0 : -1))

Definition at line 20 of file zeroInterp.c.

Function Documentation

zeroInterp()

double zeroInterp	(	double(*	fn )(),
		double	value,
		double	x_i,
		double	x_f,
		double	dx,
		double	_zero )

Finds the zero of a function within a specified interval using successive interpolation.

This function attempts to find a value x such that fn(x) is approximately equal to the given value. It uses an iterative interpolation method to locate the zero within the interval [x_i, x_f].

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 value of the independent variable (start of the interval).
x_f	Final value of the independent variable (end of the interval).
dx	Step size to use when searching the interval.
_zero	Acceptable tolerance for the zero, determining when to stop the interpolation.

Returns

The x value where fn(x) is approximately equal to value, or x_f + dx if no zero is found within the interval.

Definition at line 37 of file zeroInterp.c.

{
  double xa, xb, xm, x_b;
  double fa, fb, fm;
  double f_abs, f_bdd;
  long s_fa, s_fb, s_fm;
 
  if (dx > (x_f - x_i))
    dx = (x_f - x_i) / 2;
 
  xa = x_i;
  xb = xa + dx;
 
  if (xb > x_f)
    xb = x_f;
  if (xa == xb)
    xa = xb - dx;
 
  fa = (*fn)(xa)-value;
  s_fa = sign(fa);
 
  while (xb <= x_f) {
    fb = (*fn)(xb)-value;
    s_fb = sign(fb);
    if (s_fb == s_fa) {
      fa = fb;
      xa = xb;
      s_fa = s_fb;
      xb = xb + dx;
    } else {
      /* function has passed through zero    */
      /* so interpolate to find zero, repeat */
      f_bdd = 1000 * fabs(fa);
      xm = xa - fa * (xa - xb) / (fa - fb);
      fm = (*fn)(xm)-value;
      s_fm = sign(fm);
      x_b = xb;
      do {
        if (s_fm == 0)
          return (xm);
        else if (s_fm != s_fa) {
          xb = xm;
          fb = fm;
          s_fb = s_fm;
        } else {
          xa = xm;
          fa = fm;
          s_fa = s_fm;
        }
        xm = xa - fa * (xa - xb) / (fa - fb);
        fm = (*fn)(xm)-value;
        s_fm = sign(fm);
        f_abs = fabs(fm);
      } while (f_abs > _zero && f_abs < f_bdd);
      if (f_abs < _zero)
        return (xm);
      /* Function had a tan(x)-like singularity, which
             * looked like a zero.  So step beyond singularity
             * and continue search.
             */
      return (zeroInterp(fn, value, x_b, x_f, dx, _zero));
    }
  }
  return (x_f + dx);
}