Changeset 983 for trunk/test_more
 Timestamp:
 Aug 8, 2007 2:07:48 AM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test_more/ode_err_control.cpp
r982 r983 213 213 return ok; 214 214 } 215 215 /* 216 Define 217 $latex X : \R \rightarrow \R^2$$ by 218 $latex \[ 219 \begin{array}{rcl} 220 X_0 (0) & = & 1 \\ 221 X_1 (0) & = & 0 \\ 222 X_0^{(1)} (t) & = & 2 \alpha t X_0 (t) \\ 223 X_1^{(1)} (t) & = & 1 / X_0 (t) 224 \end{array} 225 \] $$ 226 It follows that 227 $latex \[ 228 \begin{array}{rcl} 229 X_0 (t) & = & \exp ( \alpha t^2 ) \\ 230 X_1 (t) & = & \int_0^t \exp(  \alpha s^2 ) {\bf d} s \\ 231 & = & 232 \frac{1}{ \sqrt{\alpha} \int_0^{\sqrt{\alpha} t} \exp(  r^2 ) {\bf d} r 233 \\ 234 & = & \frac{\sqrt{\pi}}{ 2 \sqrt{\alpha} {\rm erf} ( \sqrt{\alpha} t ) 235 \end{array} 236 \] $$ 237 If $latex X_0 (t) < 0$$, 238 we return $code nan$$ in order to inform 239 $code OdeErrControl$$ that its is taking to large a step. 240 241 */ 242 243 # include <cppad/rosen_34.hpp> // CppAD::Rosen34 244 # include <cppad/cppad.h> 245 246 namespace { 247 //  248 class Fun_three { 249 private: 250 const double alpha_; 251 bool was_negative_; 252 public: 253 // constructor 254 Fun_three(double alpha) : alpha_(alpha), was_negative_(false) 255 { } 256 257 // set f = x'(t) 258 void Ode( 259 const double &t, 260 const CppAD::vector<double> &x, 261 CppAD::vector<double> &f) 262 { f[0] = 2. * alpha_ * t * x[0]; 263 f[1] = 1. / x[0]; 264 // case where ODE does not make sense 265 if( x[0] < 0.  x[1] < 0. ) 266 { was_negative_ = true; 267 f[0] = CppAD::nan(0.); 268 } 269 } 270 // set f_t = df / dt 271 void Ode_ind( 272 const double &t, 273 const CppAD::vector<double> &x, 274 CppAD::vector<double> &f_t) 275 { 276 f_t[0] = 2. * alpha_ * x[0]; 277 f_t[1] = 0.; 278 if( x[0] < 0.  x[1] < 0. ) 279 { was_negative_ = true; 280 f_t[0] = CppAD::nan(0.); 281 } 282 } 283 // set f_x = df / dx 284 void Ode_dep( 285 const double &t, 286 const CppAD::vector<double> &x, 287 CppAD::vector<double> &f_x) 288 { double x0_sq = x[0] * x[0]; 289 f_x[0 * 2 + 0] = 2. * alpha_ * t; // f0 w.r.t. x0 290 f_x[0 * 2 + 1] = 0.; // f0 w.r.t. x1 291 f_x[1 * 2 + 0] = 1./x0_sq; // f1 w.r.t. x0 292 f_x[1 * 2 + 1] = 0.; // f1 w.r.t. x1 293 if( x[0] < 0.  x[1] < 0. ) 294 { was_negative_ = true; 295 f_x[0] = CppAD::nan(0.); 296 } 297 } 298 bool was_negative(void) 299 { return was_negative_; } 300 301 }; 302 303 //  304 class Method_three { 305 public: 306 Fun_three F; 307 308 // constructor 309 Method_three(double alpha) : F(alpha) 310 { } 311 void step( 312 double ta, 313 double tb, 314 CppAD::vector<double> &xa , 315 CppAD::vector<double> &xb , 316 CppAD::vector<double> &eb ) 317 { xb = CppAD::Rosen34(F, 1, ta, tb, xa, eb); 318 } 319 size_t order(void) 320 { return 3; } 321 }; 322 } 323 324 bool OdeErrControl_three(void) 325 { bool ok = true; // initial return value 326 327 double alpha = 10.; 328 Method_three method(alpha); 329 330 CppAD::vector<double> xi(2); 331 xi[0] = 1.; 332 xi[1] = 0.; 333 334 CppAD::vector<double> eabs(2); 335 eabs[0] = 1e4; 336 eabs[1] = 1e4; 337 338 // inputs 339 double ti = 0.; 340 double tf = 1.; 341 double smin = 1e4; 342 double smax = 1.; 343 double scur = 1.; 344 double erel = 0.; 345 346 // outputs 347 CppAD::vector<double> ef(2); 348 CppAD::vector<double> xf(2); 349 CppAD::vector<double> maxabs(2); 350 size_t nstep; 351 352 353 xf = OdeErrControl(method, 354 ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep); 355 356 357 double x0 = exp( alpha * tf * tf ); 358 ok &= CppAD::NearEqual(x0, xf[0], 1e4, 1e4); 359 ok &= CppAD::NearEqual(0., ef[0], 1e4, 1e4); 360 361 double root_pi = sqrt( 4. * atan(1.)); 362 double root_alpha = sqrt( alpha ); 363 double x1 = CppAD::erf(alpha * tf) * root_pi / (2 * root_alpha); 364 ok &= CppAD::NearEqual(x1, xf[1], 1e4, 1e4); 365 ok &= CppAD::NearEqual(0., ef[1], 1e4, 1e4); 366 367 ok &= method.F.was_negative(); 368 369 return ok; 370 } 371 372 // ========================================================================== 216 373 bool OdeErrControl(void) 217 374 { bool ok = true; 218 375 ok &= OdeErrControl_one(); 219 376 ok &= OdeErrControl_two(); 377 ok &= OdeErrControl_three(); 220 378 return ok; 221 379 }
Note: See TracChangeset
for help on using the changeset viewer.