cloudy  trunk
iter_track.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 
4 #include "cddefines.h"
5 #include "thirdparty.h"
6 #include "iter_track.h"
7 
9 // The class iter_track and routine Amsterdam_Method were derived from this //
10 // original source: http://www.mymathlib.webtrellis.net/roots/amsterdam.html //
11 // //
12 // The original code was heavily modified by: Peter van Hoof, ROB //
13 // //
14 // The Amsterdam method is more commonly known as the: //
15 // van Wijngaarden-Dekker-Brent method //
17 
18 double iter_track::next_val()
19 {
20  // If the function at the left endpoint is positive, and the function //
21  // at the right endpoint is negative. Iterate reducing the length //
22  // of the interval by either bisection or quadratic inverse //
23  // interpolation. Note that the function at the left endpoint //
24  // remains nonnegative and the function at the right endpoint remains //
25  // nonpositive. //
26 
27  if( p_y(p_a) > 0.0 )
28  {
29  // Check that we are converging or that we have converged near //
30  // the left endpoint of the inverval. If it appears that the //
31  // interval is not decreasing fast enough, use bisection. //
32  if( (p_x(p_b)-p_x(p_a)) < p_tol )
33  {
34  if( p_y(p_b) > 0 )
35  p_a = p_b;
36  else
37  p_set_root(p_x(p_b));
38  return p_midpoint();
39  }
40 
41  // Check that we are converging or that we have converged near //
42  // the right endpoint of the inverval. If it appears that the //
43  // interval is not decreasing fast enough, use bisection. //
44  if( (p_x(p_c)-p_x(p_b)) < p_tol )
45  {
46  if( p_y(p_b) < 0 )
47  p_c = p_b;
48  else
49  p_set_root(p_x(p_b));
50  return p_midpoint();
51  }
52 
53  // If quadratic inverse interpolation is feasible, try it. //
54 
55  if( ( p_y(p_a) > p_y(p_b) ) && ( p_y(p_b) > p_y(p_c) ) )
56  {
57  double delta = p_denominator(p_y(p_a),p_y(p_b),p_y(p_c));
58  if( delta != 0.0 )
59  {
60  double dab = p_x(p_a)-p_x(p_b);
61  double dcb = p_x(p_c)-p_x(p_b);
62  delta = safe_div( p_numerator(dab,dcb,p_y(p_a),p_y(p_b),p_y(p_c)), delta );
63 
64  // Will the new estimate of the root be within the //
65  // interval? If yes, use it and update interval. //
66  // If no, use the bisection method. //
67 
68  if( delta > dab && delta < dcb )
69  {
70  if( p_y(p_b) > 0.0 )
71  p_a = p_b;
72  else if( p_y(p_b) < 0.0 )
73  p_c = p_b;
74  else
75  p_set_root(p_x(p_b));
76  return p_x(p_b) + delta;
77  }
78  }
79  }
80 
81  // If not, use the bisection method. //
82 
83  if( p_y(p_b) > 0.0 )
84  p_a = p_b;
85  else
86  p_c = p_b;
87  return p_midpoint();
88  }
89  else
90  {
91  // If the function at the left endpoint is negative, and the function //
92  // at the right endpoint is positive. Iterate reducing the length //
93  // of the interval by either bisection or quadratic inverse //
94  // interpolation. Note that the function at the left endpoint //
95  // remains nonpositive and the function at the right endpoint remains //
96  // nonnegative. //
97 
98  if( (p_x(p_b)-p_x(p_a)) < p_tol )
99  {
100  if( p_y(p_b) < 0 )
101  p_a = p_b;
102  else
103  p_set_root(p_x(p_b));
104  return p_midpoint();
105  }
106 
107  if( (p_x(p_c)-p_x(p_b)) < p_tol )
108  {
109  if( p_y(p_b) > 0 )
110  p_c = p_b;
111  else
112  p_set_root(p_x(p_b));
113  return p_midpoint();
114  }
115 
116  if( ( p_y(p_a) < p_y(p_b) ) && ( p_y(p_b) < p_y(p_c) ) )
117  {
118  double delta = p_denominator(p_y(p_a),p_y(p_b),p_y(p_c));
119  if( delta != 0.0 )
120  {
121  double dab = p_x(p_a)-p_x(p_b);
122  double dcb = p_x(p_c)-p_x(p_b);
123  delta = safe_div( p_numerator(dab,dcb,p_y(p_a),p_y(p_b),p_y(p_c)), delta );
124  if( delta > dab && delta < dcb )
125  {
126  if( p_y(p_b) < 0.0 )
127  p_a = p_b;
128  else if( p_y(p_b) > 0.0 )
129  p_c = p_b;
130  else
131  p_set_root(p_x(p_b));
132  return p_x(p_b) + delta;
133  }
134  }
135  }
136 
137  if( p_y(p_b) < 0.0 )
138  p_a = p_b;
139  else
140  p_c = p_b;
141  return p_midpoint();
142  }
143 }
144 
145 double iter_track::deriv(int n, double& sigma) const
146 {
147  n = min( n, p_history.size() );
148  ASSERT( n >= 2 );
149  double *x = new double[n];
150  double *y = new double[n];
151  for( int i=0; i < n; ++i )
152  {
153  x[i] = p_x(p_history.size() - n + i);
154  y[i] = p_y(p_history.size() - n + i);
155  }
156  double a,b,siga,sigb;
157  linfit( n, x, y, a, siga, b, sigb );
158  delete[] y;
159  delete[] x;
160  sigma = sigb;
161  return b;
162 }
163 
164 double iter_track::zero_fit(int n, double& sigma) const
165 {
166  n = min( n, p_history.size() );
167  ASSERT( n >= 2 );
168  double *x = new double[n];
169  double *y = new double[n];
170  for( int i=0; i < n; ++i )
171  {
172  x[i] = p_y(p_history.size() - n + i);
173  y[i] = p_x(p_history.size() - n + i);
174  }
175  double a,b,siga,sigb;
176  linfit( n, x, y, a, siga, b, sigb );
177  delete[] y;
178  delete[] x;
179  sigma = siga;
180  return a;
181 }
182 
184 // double Amsterdam_Method( double (*f)(double), double a, double fa, //
185 // double c, double fc, double tolerance, //
186 // int max_iterations, int *err) //
187 // //
188 // Description: //
189 // Estimate the root (zero) of f(x) using the Amsterdam method where //
190 // 'a' and 'c' are initial estimates which bracket the root i.e. either //
191 // f(a) > 0 and f(c) < 0 or f(a) < 0 and f(c) > 0. The iteration //
192 // terminates when the zero is constrained to be within an interval of //
193 // length < 'tolerance', in which case the value returned is the best //
194 // estimate that interval. //
195 // //
196 // The Amsterdam method is an extension of Mueller's successive bisection //
197 // and inverse quadratic interpolation. Later extended by Van Wijnaarden,//
198 // Dekker and still later by Brent. Initially, the method uses the two //
199 // bracketing endpoints and the midpoint to estimate an inverse quadratic //
200 // to interpolate for the root. The interval is successively reduced by //
201 // keeping endpoints which bracket the root and an interior point used //
202 // to estimate a quadratic. If the interior point becomes too close to //
203 // an endpoint and the function has the same sign at both points, the //
204 // interior point is chosen by the bisection method. If inverse //
205 // quadratic interpolation is not feasible, the new interior point is //
206 // chosen by bisection, and if inverse quadratic interpolation results //
207 // in a point exterior to the bracketing interval, the new interior point //
208 // is again chosen by bisection. //
209 // //
210 // Arguments: //
211 // double *f Pointer to function of a single variable of type //
212 // double. //
213 // double a Initial estimate. //
214 // double fa function value at a //
215 // double c Initial estimate. //
216 // double fc function value at c //
217 // double tolerance Desired accuracy of the zero. //
218 // int max_iterations The maximum allowable number of iterations. //
219 // int *err 0 if successful, -1 if not, i.e. if f(a)*f(c) > 0, //
220 // -2 if the number of iterations > max_iterations. //
221 // //
222 // Return Values: //
223 // A zero contained within the interval (a,c). //
224 // //
225 // Example: //
226 // { //
227 // double f(double), zero, a, fa, c, fc, tol = 1.e-6; //
228 // int err, max_iter = 20; //
229 // //
230 // (determine lower bound, a, and upper bound, c, of a zero) //
231 // fa = f(a); fc = f(c); //
232 // zero = Amsterdam_Method( f, a, fa, c, fc, tol, max_iter, &err); //
233 // ... //
234 // } //
235 // double f(double x) { define f } //
237 
238 double Amsterdam_Method( double (*f)(double), double a, double fa, double c, double fc,
239  double tol, int max_iter, int& err )
240 {
241  iter_track track;
242 
243  double result;
244  set_NaN( result );
245 
246  track.set_tol(tol);
247 
248  // If the initial estimates do not bracket a root, set the err flag. //
249  if( ( err = track.init_bracket( a, fa, c, fc ) ) < 0 )
250  return result;
251 
252  double b = 0.5*(a + c);
253  for( int i = 0; i < max_iter && !track.lgConverged(); i++ )
254  {
255  track.add( b, (*f)(b) );
256  b = track.next_val();
257  }
258 
259  if( track.lgConverged() )
260  {
261  err = 0;
262  result = track.root();
263  }
264  else
265  {
266  err = -2;
267  }
268  return result;
269 }
iter_track::p_denominator
double p_denominator(double fa, double fb, double fc)
Definition: iter_track.h:62
iter_track::p_history
vector< pair< double, double > > p_history
Definition: iter_track.h:19
iter_track::p_c
int p_c
Definition: iter_track.h:24
thirdparty.h
Amsterdam_Method
double Amsterdam_Method(double(*f)(double), double a, double fa, double c, double fc, double tol, int max_iter, int &err)
Definition: iter_track.cpp:238
iter_track::p_x
double p_x(int ip) const
Definition: iter_track.h:46
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iter_track::add
void add(double x, double fx)
Definition: iter_track.h:120
iter_track::p_set_root
void p_set_root(double x)
Definition: iter_track.h:41
iter_track::root
double root() const
Definition: iter_track.h:100
iter_track::p_midpoint
double p_midpoint() const
Definition: iter_track.h:54
iter_track::lgConverged
bool lgConverged()
Definition: iter_track.h:89
iter_track::zero_fit
double zero_fit() const
Definition: iter_track.h:165
cddefines.h
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:961
iter_track::p_b
int p_b
Definition: iter_track.h:23
iter_track::deriv
double deriv() const
Definition: iter_track.h:148
set_NaN
void set_NaN(sys_float &x)
Definition: cpu.cpp:682
iter_track::init_bracket
int init_bracket(double x1, double fx1, double x2, double fx2)
Definition: iter_track.h:104
iter_track::set_tol
void set_tol(double tol)
Definition: iter_track.h:81
iter_track
Definition: iter_track.h:17
iter_track::p_numerator
double p_numerator(double dab, double dcb, double fa, double fb, double fc)
Definition: iter_track.h:58
iter_track.h
min
long min(int a, long b)
Definition: cddefines.h:723
linfit
bool linfit(long n, const double xorg[], const double yorg[], double &a, double &siga, double &b, double &sigb)
Definition: thirdparty.cpp:46
iter_track::p_a
int p_a
Definition: iter_track.h:22
iter_track::next_val
double next_val()
Definition: iter_track.cpp:18
iter_track::p_y
double p_y(int ip) const
Definition: iter_track.h:50
iter_track::p_tol
double p_tol
Definition: iter_track.h:21
Generated by   doxygen 1.8.17