cloudy trunk
Loading...
Searching...
No Matches
TestIterTrack.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#include "cdstd.h"
4#include <UnitTest++.h>
5#include "cddefines.h"
6#include "iter_track.h"
7
8namespace {
9 // test a function that would have gone into a limit cycle
10 TEST(IterTrackBasicFloat)
11 {
13 sys_float x = 1.f;
14 sys_float xnew = 2.f;
15 while( abs(x-xnew) > 2.f*FLT_EPSILON*abs(x) )
16 {
17 x = xnew;
18 // derivative at root is -1
19 xnew = track.next_val( x, 2.f/x );
20 }
21 CHECK( fp_equal( xnew, (sys_float)sqrt(2.) ) );
22 // now clear the tracker and try again
23 track.clear();
24 x = 1.f;
25 xnew = 2.f;
26 while( abs(x-xnew) > 2.f*FLT_EPSILON*abs(x) )
27 {
28 x = xnew;
29 // derivative at root is -1
30 xnew = track.next_val( x, 3.f/x );
31 }
32 CHECK( fp_equal( xnew, (sys_float)sqrt(3.) ) );
33 }
34
35 // now test the same for doubles
36 TEST(IterTrackBasicDouble)
37 {
39 double x = 1.;
40 double xnew = 2.;
41 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
42 {
43 x = xnew;
44 // derivative at root is -1
45 xnew = track.next_val( x, 2./x );
46 }
47 CHECK( fp_equal( xnew, sqrt(2.) ) );
48 }
49
50 // test a function that would have diverged
51 TEST(IterTrackBasicUnstable)
52 {
54 double x = 1.;
55 double xnew = 2.;
56 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
57 {
58 x = xnew;
59 // derivative at root is -5
60 xnew = track.next_val( x, 1./x - 2.*x );
61 }
62 // possible solutions are +/-sqrt(1/3), either one may be found
63 CHECK( fp_equal( abs(xnew), 1./sqrt(3.) ) );
64 }
65
66 // test a function that would have converged anyway
67 TEST(IterTrackBasicStableNeg)
68 {
70 double x = 1.;
71 double xnew = 2.;
72 // this would have diverged without iter_tracking
73 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
74 {
75 x = xnew;
76 // derivative at root is -1/3
77 xnew = track.next_val( x, 1./x + x/3. );
78 }
79 CHECK( fp_equal( xnew, sqrt(1.5) ) );
80 }
81
82 TEST(IterTrackBasicStablePos)
83 {
85 double x = 1.;
86 double xnew = 2.;
87 // this would have diverged without iter_tracking
88 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
89 {
90 x = xnew;
91 // derivative at root is 1/3
92 xnew = track.next_val( x, 1./x + 2.*x/3. );
93 }
94 CHECK( fp_equal( xnew, sqrt(3.) ) );
95 }
96
97 double testfun(double x)
98 {
99 return sin(x)-0.5;
100 }
101
102 TEST(IterTrack)
103 {
104 double x1, fx1, x2, fx2, x3, fx3;
105 iter_track track;
106 x1 = 0.;
107 fx1 = testfun(x1);
108 x3 = 1.5;
109 fx3 = testfun(x3);
110 double tol = 1.e-12;
111 track.set_tol(tol);
112 CHECK_EQUAL( 0, track.init_bracket(x1,fx1,x3,fx3) );
113 CHECK( fp_equal( track.bracket_width(), abs(x3-x1) ) );
114 CHECK( !track.lgConverged() );
115 x2 = 0.5*(x1+x3);
116 fx2 = testfun(x2);
117 track.add( x2, fx2 );
118 double xnew = track.next_val(0.01);
119 CHECK( fp_equal_tol( abs(xnew/x2-1.), 0.01, 1.e-12 ) );
120 const int navg = 5;
121 vector<double> xvals( navg ); // keep track of the last navg x-values
122 for( int i=0; i < 100 && !track.lgConverged(); ++i )
123 {
124 x2 = track.next_val();
125 fx2 = testfun(x2);
126 track.add( x2, fx2 );
127 // use xvals as circular buffer
128 xvals[i%navg] = x2;
129 }
130 CHECK( track.lgConverged() );
131 double exact_root = asin(0.5);
132 CHECK( fp_equal_tol( track.root(), exact_root, tol ) );
133 double sigma;
134 double val = track.deriv( navg, sigma );
135 double delta_lo = *min_element( xvals.begin(), xvals.end() ) - exact_root;
136 double delta_hi = *max_element( xvals.begin(), xvals.end() ) - exact_root;
137 CHECK( delta_lo < 0. );
138 CHECK( delta_hi > 0. );
139 // the exact derivative at the root is sqrt(3)/2 = 0.8660254...
140 // the exact 2nd derivative at the root is -1/2
141 double err_lo = -0.5*delta_lo;
142 double err_hi = -0.5*delta_hi;
143 CHECK( fp_bound( sqrt(3.)/2.+err_hi, val, sqrt(3.)/2.+err_lo ) );
144 // the tangent at the exact root is given by asin(0.5) + sqrt(3)/2*(x-x0)
145 // if we subtract that from the Taylor expansion of testfun we get:
146 // residual = -1/4*(x-x0)^2 + O((x-x0)^3), hence sigma should be less
147 // than the maximum absolute value of -1/4*(x-x0)^2 (the actual fit
148 // should run slightly closer to the maximum deviant value than the tangent).
149 CHECK( sigma < max( pow2(err_lo), pow2(err_hi) ) );
150 // ask for more points than are available to see if that is handled correctly
151 val = track.deriv( 200 );
152 double val2 = track.deriv();
153 CHECK( fp_equal( val, val2 ) );
154 val = track.deriv( 200, sigma );
155 double sigma2;
156 val = track.deriv( sigma2 );
157 CHECK( fp_equal( sigma, sigma2 ) );
158
159 // now do the same thing for the zero_fit() methods...
160 val = track.zero_fit( navg, sigma );
161 // the exact root is asin(0.5) = 0.52359877...
162 CHECK( fp_equal_tol( exact_root, val, 2.*sigma ) );
163 val = track.zero_fit( 200 );
164 val2 = track.zero_fit();
165 CHECK( fp_equal( val, val2 ) );
166 val = track.zero_fit( 200, sigma );
167 val = track.zero_fit( sigma2 );
168 CHECK( fp_equal( sigma, sigma2 ) );
169 }
170
171 // this is the short version of the above...
172 TEST(AmsterdamMethod)
173 {
174 double x1, fx1, x2, fx2;
175 x1 = 0.;
176 fx1 = testfun(x1);
177 x2 = 1.5;
178 fx2 = testfun(x2);
179 double tol = 1.e-12;
180 int err = -1;
181 double x = Amsterdam_Method( testfun, x1, fx1, x2, fx2, tol, 1000, err );
182 CHECK_EQUAL( 0, err );
183 CHECK( fp_equal_tol( x, asin(0.5), tol ) );
184 }
185
186 // now test some unstable functions
187 double testfun2(double x)
188 {
189 return exp(x)-3.;
190 }
191
192 // the derivative at the root is 3
193 TEST(AmsterdamMethod2)
194 {
195 double x1, fx1, x2, fx2;
196 x1 = 0.;
197 fx1 = testfun2(x1);
198 x2 = 3.;
199 fx2 = testfun2(x2);
200 double tol = 1.e-12;
201 int err = -1;
202 double x = Amsterdam_Method( testfun2, x1, fx1, x2, fx2, tol, 1000, err );
203 CHECK_EQUAL( 0, err );
204 CHECK( fp_equal_tol( x, log(3.), tol ) );
205 }
206
207 double testfun3(double x)
208 {
209 return 1./x - 2.*x;
210 }
211
212 // the derivative at the root is -4
213 TEST(AmsterdamMethod3)
214 {
215 double x1, fx1, x2, fx2;
216 x1 = 0.1;
217 fx1 = testfun3(x1);
218 x2 = 3.;
219 fx2 = testfun3(x2);
220 double tol = 1.e-12;
221 int err = -1;
222 double x = Amsterdam_Method( testfun3, x1, fx1, x2, fx2, tol, 1000, err );
223 CHECK_EQUAL( 0, err );
224 CHECK( fp_equal_tol( x, sqrt(0.5), tol ) );
225 }
226}
static double x2[63]
static double x1[83]
float sys_float
Definition cddefines.h:106
T pow2(T a)
Definition cddefines.h:931
bool fp_bound(sys_float lo, sys_float x, sys_float hi, int n=3)
Definition cddefines.h:877
long max(int a, long b)
Definition cddefines.h:775
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition cddefines.h:854
T next_val(T current, T next_est)
Definition iter_track.h:250
void add(double x, double fx)
Definition iter_track.h:120
double bracket_width() const
Definition iter_track.h:85
bool lgConverged()
Definition iter_track.h:89
int init_bracket(double x1, double fx1, double x2, double fx2)
Definition iter_track.h:104
double root() const
Definition iter_track.h:100
double zero_fit(int n, double &sigma) const
double next_val()
void set_tol(double tol)
Definition iter_track.h:81
double deriv(int n, double &sigma) const
double Amsterdam_Method(double(*f)(double), double a, double fa, double c, double fc, double tol, int max_iter, int &err)