cloudy trunk
Loading...
Searching...
No Matches
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
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
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
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
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
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
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
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
145double 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
164double 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
238double 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}
#define ASSERT(exp)
Definition cddefines.h:578
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition cddefines.h:961
long min(int a, long b)
Definition cddefines.h:723
double p_numerator(double dab, double dcb, double fa, double fb, double fc)
Definition iter_track.h:58
double deriv() const
Definition iter_track.h:148
void add(double x, double fx)
Definition iter_track.h:120
double p_denominator(double fa, double fb, double fc)
Definition iter_track.h:62
double p_tol
Definition iter_track.h:21
bool lgConverged()
Definition iter_track.h:89
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
double root() const
Definition iter_track.h:100
double next_val()
double p_y(int ip) const
Definition iter_track.h:50
vector< pair< double, double > > p_history
Definition iter_track.h:19
double p_x(int ip) const
Definition iter_track.h:46
void set_tol(double tol)
Definition iter_track.h:81
double p_midpoint() const
Definition iter_track.h:54
void p_set_root(double x)
Definition iter_track.h:41
void set_NaN(sys_float &x)
Definition cpu.cpp:682
double Amsterdam_Method(double(*f)(double), double a, double fa, double c, double fc, double tol, int max_iter, int &err)
bool linfit(long n, const double xorg[], const double yorg[], double &a, double &siga, double &b, double &sigb)