cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iter_track.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 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 "iter_track.h"
6 #include "thirdparty.h"
7 
9 // The class iter_track and routine Amsterdam_Method were derived from this //
10 // original source: http://www.mymathlib.com/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 
20 {
21  realnum diff = next_est-current;
22  realnum proposal = next_est;
23  if (m_lastnext != 0.0 && diff != 0.0)
24  {
25  realnum lastdiff = m_lastnext-m_lastcurr;
26  realnum numer = m_lastcurr*next_est-m_lastnext*current;
27  realnum denom = diff-lastdiff;
28  if (fabs(numer) >= 2.0*fabs(current*denom))
29  {
30  m_confidence *= 0.5;
31  }
32  else
33  {
34  proposal = numer/denom;
35  //fprintf(ioQQQ,"!!!! %g %g %g %g\n",proposal,current,next_est,diff);
36  //fflush(ioQQQ);
37  if ((proposal-current)*(proposal-next_est) > 0.2*diff*diff)
38  {
39  proposal = next_est;
40  m_confidence *= 0.5;
41  }
42  else
43  {
44  m_confidence *= 2.0;
45  }
46  }
47  m_confidence = min(1.0,max(0.25,m_confidence));
48 
49  proposal = m_confidence*proposal+(1.0-m_confidence)*current;
50  }
51  m_lastnext = next_est;
52  m_lastcurr = current;
53  return proposal;
54 }
55 
56 
58 {
59  // If the function at the left endpoint is positive, and the function //
60  // at the right endpoint is negative. Iterate reducing the length //
61  // of the interval by either bisection or quadratic inverse //
62  // interpolation. Note that the function at the left endpoint //
63  // remains nonnegative and the function at the right endpoint remains //
64  // nonpositive. //
65 
66  if( p_y(p_a) > 0.0 )
67  {
68  // Check that we are converging or that we have converged near //
69  // the left endpoint of the inverval. If it appears that the //
70  // interval is not decreasing fast enough, use bisection. //
71  if( (p_x(p_b)-p_x(p_a)) < p_tol )
72  {
73  if( p_y(p_b) > 0 )
74  p_a = p_b;
75  else
76  p_set_root(p_x(p_b));
77  return p_midpoint();
78  }
79 
80  // Check that we are converging or that we have converged near //
81  // the right endpoint of the inverval. If it appears that the //
82  // interval is not decreasing fast enough, use bisection. //
83  if( (p_x(p_c)-p_x(p_b)) < p_tol )
84  {
85  if( p_y(p_b) < 0 )
86  p_c = p_b;
87  else
88  p_set_root(p_x(p_b));
89  return p_midpoint();
90  }
91 
92  // If quadratic inverse interpolation is feasible, try it. //
93 
94  if( ( p_y(p_a) > p_y(p_b) ) && ( p_y(p_b) > p_y(p_c) ) )
95  {
96  double delta = p_denominator(p_y(p_a),p_y(p_b),p_y(p_c));
97  if( delta != 0.0 )
98  {
99  double dab = p_x(p_a)-p_x(p_b);
100  double dcb = p_x(p_c)-p_x(p_b);
101  delta = safe_div( p_numerator(dab,dcb,p_y(p_a),p_y(p_b),p_y(p_c)), delta );
102 
103  // Will the new estimate of the root be within the //
104  // interval? If yes, use it and update interval. //
105  // If no, use the bisection method. //
106 
107  if( delta > dab && delta < dcb )
108  {
109  if( p_y(p_b) > 0.0 )
110  p_a = p_b;
111  else if( p_y(p_b) < 0.0 )
112  p_c = p_b;
113  else
114  p_set_root(p_x(p_b));
115  return p_x(p_b) + delta;
116  }
117  }
118  }
119 
120  // If not, use the bisection method. //
121 
122  if( p_y(p_b) > 0.0 )
123  p_a = p_b;
124  else
125  p_c = p_b;
126  return p_midpoint();
127  }
128  else
129  {
130  // If the function at the left endpoint is negative, and the function //
131  // at the right endpoint is positive. Iterate reducing the length //
132  // of the interval by either bisection or quadratic inverse //
133  // interpolation. Note that the function at the left endpoint //
134  // remains nonpositive and the function at the right endpoint remains //
135  // nonnegative. //
136 
137  if( (p_x(p_b)-p_x(p_a)) < p_tol )
138  {
139  if( p_y(p_b) < 0 )
140  p_a = p_b;
141  else
142  p_set_root(p_x(p_b));
143  return p_midpoint();
144  }
145 
146  if( (p_x(p_c)-p_x(p_b)) < p_tol )
147  {
148  if( p_y(p_b) > 0 )
149  p_c = p_b;
150  else
151  p_set_root(p_x(p_b));
152  return p_midpoint();
153  }
154 
155  if( ( p_y(p_a) < p_y(p_b) ) && ( p_y(p_b) < p_y(p_c) ) )
156  {
157  double delta = p_denominator(p_y(p_a),p_y(p_b),p_y(p_c));
158  if( delta != 0.0 )
159  {
160  double dab = p_x(p_a)-p_x(p_b);
161  double dcb = p_x(p_c)-p_x(p_b);
162  delta = safe_div( p_numerator(dab,dcb,p_y(p_a),p_y(p_b),p_y(p_c)), delta );
163  if( delta > dab && delta < dcb )
164  {
165  if( p_y(p_b) < 0.0 )
166  p_a = p_b;
167  else if( p_y(p_b) > 0.0 )
168  p_c = p_b;
169  else
170  p_set_root(p_x(p_b));
171  return p_x(p_b) + delta;
172  }
173  }
174  }
175 
176  if( p_y(p_b) < 0.0 )
177  p_a = p_b;
178  else
179  p_c = p_b;
180  return p_midpoint();
181  }
182 }
183 
184 double iter_track::deriv(int n, double& sigma) const
185 {
186  n = min( n, p_history.size() );
187  ASSERT( n >= 2 );
188  double *x = new double[n];
189  double *y = new double[n];
190  for( int i=0; i < n; ++i )
191  {
192  x[i] = p_x(p_history.size() - n + i);
193  y[i] = p_y(p_history.size() - n + i);
194  }
195  double a,b,siga,sigb;
196  linfit( n, x, y, a, siga, b, sigb );
197  delete[] y;
198  delete[] x;
199  sigma = sigb;
200  return b;
201 }
202 
203 double iter_track::zero_fit(int n, double& sigma) const
204 {
205  n = min( n, p_history.size() );
206  ASSERT( n >= 2 );
207  double *x = new double[n];
208  double *y = new double[n];
209  for( int i=0; i < n; ++i )
210  {
211  x[i] = p_y(p_history.size() - n + i);
212  y[i] = p_x(p_history.size() - n + i);
213  }
214  double a,b,siga,sigb;
215  linfit( n, x, y, a, siga, b, sigb );
216  delete[] y;
217  delete[] x;
218  sigma = siga;
219  return a;
220 }
221 
223 // double Amsterdam_Method( double (*f)(double), double a, double fa, //
224 // double c, double fc, double tolerance, //
225 // int max_iterations, int *err) //
226 // //
227 // Description: //
228 // Estimate the root (zero) of f(x) using the Amsterdam method where //
229 // 'a' and 'c' are initial estimates which bracket the root i.e. either //
230 // f(a) > 0 and f(c) < 0 or f(a) < 0 and f(c) > 0. The iteration //
231 // terminates when the zero is constrained to be within an interval of //
232 // length < 'tolerance', in which case the value returned is the best //
233 // estimate that interval. //
234 // //
235 // The Amsterdam method is an extension of Mueller's successive bisection //
236 // and inverse quadratic interpolation. Later extended by Van Wijnaarden,//
237 // Dekker and still later by Brent. Initially, the method uses the two //
238 // bracketing endpoints and the midpoint to estimate an inverse quadratic //
239 // to interpolate for the root. The interval is successively reduced by //
240 // keeping endpoints which bracket the root and an interior point used //
241 // to estimate a quadratic. If the interior point becomes too close to //
242 // an endpoint and the function has the same sign at both points, the //
243 // interior point is chosen by the bisection method. If inverse //
244 // quadratic interpolation is not feasible, the new interior point is //
245 // chosen by bisection, and if inverse quadratic interpolation results //
246 // in a point exterior to the bracketing interval, the new interior point //
247 // is again chosen by bisection. //
248 // //
249 // Arguments: //
250 // double *f Pointer to function of a single variable of type //
251 // double. //
252 // double a Initial estimate. //
253 // double fa function value at a //
254 // double c Initial estimate. //
255 // double fc function value at c //
256 // double tolerance Desired accuracy of the zero. //
257 // int max_iterations The maximum allowable number of iterations. //
258 // int *err 0 if successful, -1 if not, i.e. if f(a)*f(c) > 0, //
259 // -2 if the number of iterations > max_iterations. //
260 // //
261 // Return Values: //
262 // A zero contained within the interval (a,c). //
263 // //
264 // Example: //
265 // { //
266 // double f(double), zero, a, fa, c, fc, tol = 1.e-6; //
267 // int err, max_iter = 20; //
268 // //
269 // (determine lower bound, a, and upper bound, c, of a zero) //
270 // fa = f(a); fc = f(c); //
271 // zero = Amsterdam_Method( f, a, fa, c, fc, tol, max_iter, &err); //
272 // ... //
273 // } //
274 // double f(double x) { define f } //
276 
277 double Amsterdam_Method( double (*f)(double), double a, double fa, double c, double fc,
278  double tol, int max_iter, int& err )
279 {
280  iter_track track;
281 
282  double result;
283  set_NaN( result );
284 
285  track.set_tol(tol);
286 
287  // If the initial estimates do not bracket a root, set the err flag. //
288  if( ( err = track.init_bracket( a, fa, c, fc ) ) < 0 )
289  return result;
290 
291  double b = 0.5*(a + c);
292  for( int i = 0; i < max_iter && !track.lgConverged(); i++ )
293  {
294  track.add( b, (*f)(b) );
295  b = track.next_val();
296  }
297 
298  if( track.lgConverged() )
299  {
300  err = 0;
301  result = track.root();
302  }
303  else
304  {
305  err = -2;
306  }
307  return result;
308 }
double p_midpoint() const
Definition: iter_track.h:54
realnum m_lastnext
Definition: iter_track.h:270
realnum m_lastcurr
Definition: iter_track.h:271
void set_NaN(sys_float &x)
Definition: cpu.cpp:906
double next_val()
Definition: iter_track.cpp:57
void add(double x, double fx)
Definition: iter_track.h:120
double p_x(int ip) const
Definition: iter_track.h:46
double p_denominator(double fa, double fb, double fc)
Definition: iter_track.h:62
realnum m_confidence
Definition: iter_track.h:269
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:277
double p_y(int ip) const
Definition: iter_track.h:50
realnum next_val(realnum current, realnum next_est)
Definition: iter_track.cpp:19
vector< pair< double, double > > p_history
Definition: iter_track.h:19
int init_bracket(double x1, double fx1, double x2, double fx2)
Definition: iter_track.h:104
double zero_fit() const
Definition: iter_track.h:165
bool lgConverged()
Definition: iter_track.h:89
float realnum
Definition: cddefines.h:124
long max(int a, long b)
Definition: cddefines.h:817
long min(int a, long b)
Definition: cddefines.h:762
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1011
bool linfit(long n, const double xorg[], const double yorg[], double &a, double &siga, double &b, double &sigb)
Definition: thirdparty.cpp:46
#define ASSERT(exp)
Definition: cddefines.h:613
double p_tol
Definition: iter_track.h:21
void p_set_root(double x)
Definition: iter_track.h:41
double deriv() const
Definition: iter_track.h:148
void set_tol(double tol)
Definition: iter_track.h:81
double p_numerator(double dab, double dcb, double fa, double fb, double fc)
Definition: iter_track.h:58
double root() const
Definition: iter_track.h:100