00001 #include <UnitTest++.h>
00002 #include "cddefines.h"
00003 #include "iter_track.h"
00004
00005 namespace {
00006
00007 TEST(IterTrackBasicFloat)
00008 {
00009 iter_track_basic<sys_float> track;
00010 sys_float x = 1.f;
00011 sys_float xnew = 2.f;
00012 while( abs(x-xnew) > 2.f*FLT_EPSILON*abs(x) )
00013 {
00014 x = xnew;
00015
00016 xnew = track.next_val( x, 2.f/x );
00017 }
00018 CHECK( fp_equal( xnew, (sys_float)sqrt(2.) ) );
00019
00020 track.clear();
00021 x = 1.f;
00022 xnew = 2.f;
00023 while( abs(x-xnew) > 2.f*FLT_EPSILON*abs(x) )
00024 {
00025 x = xnew;
00026
00027 xnew = track.next_val( x, 3.f/x );
00028 }
00029 CHECK( fp_equal( xnew, (sys_float)sqrt(3.) ) );
00030 }
00031
00032
00033 TEST(IterTrackBasicDouble)
00034 {
00035 iter_track_basic<double> track;
00036 double x = 1.;
00037 double xnew = 2.;
00038 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
00039 {
00040 x = xnew;
00041
00042 xnew = track.next_val( x, 2./x );
00043 }
00044 CHECK( fp_equal( xnew, sqrt(2.) ) );
00045 }
00046
00047
00048 TEST(IterTrackBasicUnstable)
00049 {
00050 iter_track_basic<double> track;
00051 double x = 1.;
00052 double xnew = 2.;
00053 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
00054 {
00055 x = xnew;
00056
00057 xnew = track.next_val( x, 1./x - 2.*x );
00058 }
00059
00060 CHECK( fp_equal( abs(xnew), 1./sqrt(3.) ) );
00061 }
00062
00063
00064 TEST(IterTrackBasicStableNeg)
00065 {
00066 iter_track_basic<double> track;
00067 double x = 1.;
00068 double xnew = 2.;
00069
00070 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
00071 {
00072 x = xnew;
00073
00074 xnew = track.next_val( x, 1./x + x/3. );
00075 }
00076 CHECK( fp_equal( xnew, sqrt(1.5) ) );
00077 }
00078
00079 TEST(IterTrackBasicStablePos)
00080 {
00081 iter_track_basic<double> track;
00082 double x = 1.;
00083 double xnew = 2.;
00084
00085 while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
00086 {
00087 x = xnew;
00088
00089 xnew = track.next_val( x, 1./x + 2.*x/3. );
00090 }
00091 CHECK( fp_equal( xnew, sqrt(3.) ) );
00092 }
00093
00094 double testfun(double x)
00095 {
00096 return sin(x)-0.5;
00097 }
00098
00099 TEST(IterTrack)
00100 {
00101 double x1, fx1, x2, fx2, x3, fx3;
00102 iter_track track;
00103 x1 = 0.;
00104 fx1 = testfun(x1);
00105 x3 = 1.5;
00106 fx3 = testfun(x3);
00107 double tol = 1.e-12;
00108 track.set_tol(tol);
00109 CHECK_EQUAL( 0, track.init_bracket(x1,fx1,x3,fx3) );
00110 CHECK( fp_equal( track.bracket_width(), abs(x3-x1) ) );
00111 CHECK( !track.lgConverged() );
00112 x2 = 0.5*(x1+x3);
00113 fx2 = testfun(x2);
00114 track.add( x2, fx2 );
00115 double xnew = track.next_val(0.01);
00116 CHECK( fp_equal_tol( abs(xnew/x2-1.), 0.01, 1.e-12 ) );
00117 const int navg = 5;
00118 vector<double> xvals( navg );
00119 for( int i=0; i < 100 && !track.lgConverged(); ++i )
00120 {
00121 x2 = track.next_val();
00122 fx2 = testfun(x2);
00123 track.add( x2, fx2 );
00124
00125 xvals[i%navg] = x2;
00126 }
00127 CHECK( track.lgConverged() );
00128 double exact_root = asin(0.5);
00129 CHECK( fp_equal_tol( track.root(), exact_root, tol ) );
00130 double sigma;
00131 double val = track.deriv( navg, sigma );
00132 double delta_lo = *min_element( xvals.begin(), xvals.end() ) - exact_root;
00133 double delta_hi = *max_element( xvals.begin(), xvals.end() ) - exact_root;
00134 CHECK( delta_lo < 0. );
00135 CHECK( delta_hi > 0. );
00136
00137
00138 double err_lo = -0.5*delta_lo;
00139 double err_hi = -0.5*delta_hi;
00140 CHECK( fp_bound( sqrt(3.)/2.+err_hi, val, sqrt(3.)/2.+err_lo ) );
00141
00142
00143
00144
00145
00146 CHECK( sigma < max( pow2(err_lo), pow2(err_hi) ) );
00147
00148 val = track.deriv( 200 );
00149 double val2 = track.deriv();
00150 CHECK( fp_equal( val, val2 ) );
00151 val = track.deriv( 200, sigma );
00152 double sigma2;
00153 val = track.deriv( sigma2 );
00154 CHECK( fp_equal( sigma, sigma2 ) );
00155
00156
00157 val = track.zero_fit( navg, sigma );
00158
00159 CHECK( fp_equal_tol( exact_root, val, 2.*sigma ) );
00160 val = track.zero_fit( 200 );
00161 val2 = track.zero_fit();
00162 CHECK( fp_equal( val, val2 ) );
00163 val = track.zero_fit( 200, sigma );
00164 val = track.zero_fit( sigma2 );
00165 CHECK( fp_equal( sigma, sigma2 ) );
00166 }
00167
00168
00169 TEST(AmsterdamMethod)
00170 {
00171 double x1, fx1, x2, fx2;
00172 x1 = 0.;
00173 fx1 = testfun(x1);
00174 x2 = 1.5;
00175 fx2 = testfun(x2);
00176 double tol = 1.e-12;
00177 int err = -1;
00178 double x = Amsterdam_Method( testfun, x1, fx1, x2, fx2, tol, 1000, err );
00179 CHECK_EQUAL( 0, err );
00180 CHECK( fp_equal_tol( x, asin(0.5), tol ) );
00181 }
00182
00183
00184 double testfun2(double x)
00185 {
00186 return exp(x)-3.;
00187 }
00188
00189
00190 TEST(AmsterdamMethod2)
00191 {
00192 double x1, fx1, x2, fx2;
00193 x1 = 0.;
00194 fx1 = testfun2(x1);
00195 x2 = 3.;
00196 fx2 = testfun2(x2);
00197 double tol = 1.e-12;
00198 int err = -1;
00199 double x = Amsterdam_Method( testfun2, x1, fx1, x2, fx2, tol, 1000, err );
00200 CHECK_EQUAL( 0, err );
00201 CHECK( fp_equal_tol( x, log(3.), tol ) );
00202 }
00203
00204 double testfun3(double x)
00205 {
00206 return 1./x - 2.*x;
00207 }
00208
00209
00210 TEST(AmsterdamMethod3)
00211 {
00212 double x1, fx1, x2, fx2;
00213 x1 = 0.1;
00214 fx1 = testfun3(x1);
00215 x2 = 3.;
00216 fx2 = testfun3(x2);
00217 double tol = 1.e-12;
00218 int err = -1;
00219 double x = Amsterdam_Method( testfun3, x1, fx1, x2, fx2, tol, 1000, err );
00220 CHECK_EQUAL( 0, err );
00221 CHECK( fp_equal_tol( x, sqrt(0.5), tol ) );
00222 }
00223 }