1/*
2 [auto_generated]
3 boost/numeric/odeint/stepper/bulirsch_stoer.hpp
4
5 [begin_description]
6 Implementation of the Burlish-Stoer method. As described in
7 Ernst Hairer, Syvert Paul Norsett, Gerhard Wanner
8 Solving Ordinary Differential Equations I. Nonstiff Problems.
9 Springer Series in Comput. Mathematics, Vol. 8, Springer-Verlag 1987, Second revised edition 1993.
10 [end_description]
11
12 Copyright 2011-2013 Mario Mulansky
13 Copyright 2011-2013 Karsten Ahnert
14 Copyright 2012 Christoph Koke
15
16 Distributed under the Boost Software License, Version 1.0.
17 (See accompanying file LICENSE_1_0.txt or
18 copy at http://www.boost.org/LICENSE_1_0.txt)
19*/
20
21
22#ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED
23#define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED
24
25
26#include <iostream>
27
28#include <algorithm>
29
30#include <boost/config.hpp> // for min/max guidelines
31
32#include <boost/numeric/odeint/util/bind.hpp>
33#include <boost/numeric/odeint/util/unwrap_reference.hpp>
34
35#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
36#include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
37#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
38#include <boost/numeric/odeint/algebra/range_algebra.hpp>
39#include <boost/numeric/odeint/algebra/default_operations.hpp>
40#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
41#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
42
43#include <boost/numeric/odeint/util/state_wrapper.hpp>
44#include <boost/numeric/odeint/util/is_resizeable.hpp>
45#include <boost/numeric/odeint/util/resizer.hpp>
46#include <boost/numeric/odeint/util/unit_helper.hpp>
47#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
48
49namespace boost {
50namespace numeric {
51namespace odeint {
52
53template<
54 class State ,
55 class Value = double ,
56 class Deriv = State ,
57 class Time = Value ,
58 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
59 class Operations = typename operations_dispatcher< State >::operations_type ,
60 class Resizer = initially_resizer
61 >
62class bulirsch_stoer {
63
64public:
65
66 typedef State state_type;
67 typedef Value value_type;
68 typedef Deriv deriv_type;
69 typedef Time time_type;
70 typedef Algebra algebra_type;
71 typedef Operations operations_type;
72 typedef Resizer resizer_type;
73#ifndef DOXYGEN_SKIP
74 typedef state_wrapper< state_type > wrapped_state_type;
75 typedef state_wrapper< deriv_type > wrapped_deriv_type;
76 typedef controlled_stepper_tag stepper_category;
77
78 typedef bulirsch_stoer< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type;
79
80 typedef typename inverse_time< time_type >::type inv_time_type;
81
82 typedef std::vector< value_type > value_vector;
83 typedef std::vector< time_type > time_vector;
84 typedef std::vector< inv_time_type > inv_time_vector; //should be 1/time_type for boost.units
85 typedef std::vector< value_vector > value_matrix;
86 typedef std::vector< size_t > int_vector;
87 typedef std::vector< wrapped_state_type > state_table_type;
88#endif //DOXYGEN_SKIP
89 const static size_t m_k_max = 8;
90
91 bulirsch_stoer(
92 value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
93 value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
94 time_type max_dt = static_cast<time_type>(0))
95 : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
96 m_last_step_rejected( false ) , m_first( true ) ,
97 m_max_dt(max_dt) ,
98 m_interval_sequence( m_k_max+1 ) ,
99 m_coeff( m_k_max+1 ) ,
100 m_cost( m_k_max+1 ) ,
101 m_table( m_k_max ) ,
102 STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
103 {
104 BOOST_USING_STD_MIN();
105 BOOST_USING_STD_MAX();
106 /* initialize sequence of stage numbers and work */
107 for( unsigned short i = 0; i < m_k_max+1; i++ )
108 {
109 m_interval_sequence[i] = 2 * (i+1);
110 if( i == 0 )
111 m_cost[i] = m_interval_sequence[i];
112 else
113 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
114 m_coeff[i].resize(i);
115 for( size_t k = 0 ; k < i ; ++k )
116 {
117 const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
118 m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation
119 }
120
121 // crude estimate of optimal order
122
123 m_current_k_opt = 4;
124 /* no calculation because log10 might not exist for value_type!
125 const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 );
126 m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( m_k_max-1 ) , logfact ));
127 */
128 }
129
130 }
131
132
133 /*
134 * Version 1 : try_step( sys , x , t , dt )
135 *
136 * The overloads are needed to solve the forwarding problem
137 */
138 template< class System , class StateInOut >
139 controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt )
140 {
141 return try_step_v1( system , x , t, dt );
142 }
143
144 /**
145 * \brief Second version to solve the forwarding problem, can be used with Boost.Range as StateInOut.
146 */
147 template< class System , class StateInOut >
148 controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt )
149 {
150 return try_step_v1( system , x , t, dt );
151 }
152
153 /*
154 * Version 2 : try_step( sys , x , dxdt , t , dt )
155 *
156 * this version does not solve the forwarding problem, boost.range can not be used
157 */
158 template< class System , class StateInOut , class DerivIn >
159 controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
160 {
161 m_xnew_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_xnew< StateInOut > , detail::ref( *this ) , detail::_1 ) );
162 controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt );
163 if( res == success )
164 {
165 boost::numeric::odeint::copy( m_xnew.m_v , x );
166 }
167 return res;
168 }
169
170 /*
171 * Version 3 : try_step( sys , in , t , out , dt )
172 *
173 * this version does not solve the forwarding problem, boost.range can not be used
174 */
175 template< class System , class StateIn , class StateOut >
176 typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type
177 try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
178 {
179 typename odeint::unwrap_reference< System >::type &sys = system;
180 m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateIn > , detail::ref( *this ) , detail::_1 ) );
181 sys( in , m_dxdt.m_v , t );
182 return try_step( system , in , m_dxdt.m_v , t , out , dt );
183 }
184
185
186 /*
187 * Full version : try_step( sys , in , dxdt_in , t , out , dt )
188 *
189 * contains the actual implementation
190 */
191 template< class System , class StateIn , class DerivIn , class StateOut >
192 controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
193 {
194 if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
195 {
196 // given step size is bigger then max_dt
197 // set limit and return fail
198 dt = m_max_dt;
199 return fail;
200 }
201
202 BOOST_USING_STD_MIN();
203 BOOST_USING_STD_MAX();
204
205 static const value_type val1( 1.0 );
206
207 if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
208 {
209 reset(); // system resized -> reset
210 }
211
212 if( dt != m_dt_last )
213 {
214 reset(); // step size changed from outside -> reset
215 }
216
217 bool reject( true );
218
219 time_vector h_opt( m_k_max+1 );
220 inv_time_vector work( m_k_max+1 );
221
222 time_type new_h = dt;
223
224 /* m_current_k_opt is the estimated current optimal stage number */
225 for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
226 {
227 /* the stage counts are stored in m_interval_sequence */
228 m_midpoint.set_steps( m_interval_sequence[k] );
229 if( k == 0 )
230 {
231 m_midpoint.do_step( system , in , dxdt , t , out , dt );
232 /* the first step, nothing more to do */
233 }
234 else
235 {
236 m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt );
237 extrapolate( k , m_table , m_coeff , out );
238 // get error estimate
239 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
240 typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) );
241 const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt );
242 h_opt[k] = calc_h_opt( h: dt , error , k );
243 work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k];
244
245 if( (k == m_current_k_opt-1) || m_first )
246 { // convergence before k_opt ?
247 if( error < 1.0 )
248 {
249 //convergence
250 reject = false;
251 if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) )
252 {
253 // leave order as is (except we were in first round)
254 m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(k)+1 ) );
255 new_h = h_opt[k];
256 new_h *= static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] );
257 } else {
258 m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(k) ) );
259 new_h = h_opt[k];
260 }
261 break;
262 }
263 else if( should_reject( error , k ) && !m_first )
264 {
265 reject = true;
266 new_h = h_opt[k];
267 break;
268 }
269 }
270 if( k == m_current_k_opt )
271 { // convergence at k_opt ?
272 if( error < 1.0 )
273 {
274 //convergence
275 reject = false;
276 if( (work[k-1] < KFAC2*work[k]) )
277 {
278 m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
279 new_h = h_opt[m_current_k_opt];
280 }
281 else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected )
282 {
283 m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max-1) , static_cast<int>(m_current_k_opt)+1 );
284 new_h = h_opt[k];
285 new_h *= m_cost[m_current_k_opt]/m_cost[k];
286 } else
287 new_h = h_opt[m_current_k_opt];
288 break;
289 }
290 else if( should_reject( error , k ) )
291 {
292 reject = true;
293 new_h = h_opt[m_current_k_opt];
294 break;
295 }
296 }
297 if( k == m_current_k_opt+1 )
298 { // convergence at k_opt+1 ?
299 if( error < 1.0 )
300 { //convergence
301 reject = false;
302 if( work[k-2] < KFAC2*work[k-1] )
303 m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 );
304 if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected )
305 m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(k) );
306 new_h = h_opt[m_current_k_opt];
307 } else
308 {
309 reject = true;
310 new_h = h_opt[m_current_k_opt];
311 }
312 break;
313 }
314 }
315 }
316
317 if( !reject )
318 {
319 t += dt;
320 }
321
322 if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) )
323 {
324 // limit step size
325 if( m_max_dt != static_cast<time_type>(0) )
326 {
327 new_h = detail::min_abs(m_max_dt, new_h);
328 }
329 m_dt_last = new_h;
330 dt = new_h;
331 }
332
333 m_last_step_rejected = reject;
334 m_first = false;
335
336 if( reject )
337 return fail;
338 else
339 return success;
340 }
341
342 /** \brief Resets the internal state of the stepper */
343 void reset()
344 {
345 m_first = true;
346 m_last_step_rejected = false;
347 }
348
349
350 /* Resizer methods */
351
352 template< class StateIn >
353 void adjust_size( const StateIn &x )
354 {
355 resize_m_dxdt( x );
356 resize_m_xnew( x );
357 resize_impl( x );
358 m_midpoint.adjust_size( x );
359 }
360
361
362private:
363
364 template< class StateIn >
365 bool resize_m_dxdt( const StateIn &x )
366 {
367 return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
368 }
369
370 template< class StateIn >
371 bool resize_m_xnew( const StateIn &x )
372 {
373 return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable<state_type>::type() );
374 }
375
376 template< class StateIn >
377 bool resize_impl( const StateIn &x )
378 {
379 bool resized( false );
380 for( size_t i = 0 ; i < m_k_max ; ++i )
381 resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable<state_type>::type() );
382 resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable<state_type>::type() );
383 return resized;
384 }
385
386
387 template< class System , class StateInOut >
388 controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt )
389 {
390 typename odeint::unwrap_reference< System >::type &sys = system;
391 m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateInOut > , detail::ref( *this ) , detail::_1 ) );
392 sys( x , m_dxdt.m_v ,t );
393 return try_step( system , x , m_dxdt.m_v , t , dt );
394 }
395
396
397 template< class StateInOut >
398 void extrapolate( size_t k , state_table_type &table , const value_matrix &coeff , StateInOut &xest )
399 /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
400 uses the obtained intermediate results to extrapolate to dt->0
401 */
402 {
403 static const value_type val1 = static_cast< value_type >( 1.0 );
404 for( int j=k-1 ; j>0 ; --j )
405 {
406 m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v ,
407 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][j] , -coeff[k][j] ) );
408 }
409 m_algebra.for_each3( xest , table[0].m_v , xest ,
410 typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][0] , -coeff[k][0]) );
411 }
412
413 time_type calc_h_opt( time_type h , value_type error , size_t k ) const
414 /* calculates the optimal step size for a given error and stage number */
415 {
416 BOOST_USING_STD_MIN();
417 BOOST_USING_STD_MAX();
418 using std::pow;
419 value_type expo( 1.0/(2*k+1) );
420 value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
421 value_type fac;
422 if (error == 0.0)
423 fac=1.0/facmin;
424 else
425 {
426 fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
427 fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(facmin/STEPFAC4) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(1.0/facmin) , fac ) );
428 }
429 return h*fac;
430 }
431
432 controlled_step_result set_k_opt( size_t k , const inv_time_vector &work , const time_vector &h_opt , time_type &dt )
433 /* calculates the optimal stage number */
434 {
435 if( k == 1 )
436 {
437 m_current_k_opt = 2;
438 return success;
439 }
440 if( (work[k-1] < KFAC1*work[k]) || (k == m_k_max) )
441 { // order decrease
442 m_current_k_opt = k-1;
443 dt = h_opt[ m_current_k_opt ];
444 return success;
445 }
446 else if( (work[k] < KFAC2*work[k-1]) || m_last_step_rejected || (k == m_k_max-1) )
447 { // same order - also do this if last step got rejected
448 m_current_k_opt = k;
449 dt = h_opt[ m_current_k_opt ];
450 return success;
451 }
452 else
453 { // order increase - only if last step was not rejected
454 m_current_k_opt = k+1;
455 dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ;
456 return success;
457 }
458 }
459
460 bool in_convergence_window( size_t k ) const
461 {
462 if( (k == m_current_k_opt-1) && !m_last_step_rejected )
463 return true; // decrease stepsize only if last step was not rejected
464 return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) );
465 }
466
467 bool should_reject( value_type error , size_t k ) const
468 {
469 if( k == m_current_k_opt-1 )
470 {
471 const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] /
472 (m_interval_sequence[0]*m_interval_sequence[0]);
473 //step will fail, criterion 17.3.17 in NR
474 return ( error > d*d );
475 }
476 else if( k == m_current_k_opt )
477 {
478 const value_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0];
479 return ( error > d*d );
480 } else
481 return error > 1.0;
482 }
483
484 default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
485 modified_midpoint< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
486
487 bool m_last_step_rejected;
488 bool m_first;
489
490 time_type m_dt_last;
491 time_type m_t_last;
492 time_type m_max_dt;
493
494 size_t m_current_k_opt;
495
496 algebra_type m_algebra;
497
498 resizer_type m_dxdt_resizer;
499 resizer_type m_xnew_resizer;
500 resizer_type m_resizer;
501
502 wrapped_state_type m_xnew;
503 wrapped_state_type m_err;
504 wrapped_deriv_type m_dxdt;
505
506 int_vector m_interval_sequence; // stores the successive interval counts
507 value_matrix m_coeff;
508 int_vector m_cost; // costs for interval count
509
510 state_table_type m_table; // sequence of states for extrapolation
511
512 value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
513};
514
515
516/******** DOXYGEN ********/
517/**
518 * \class bulirsch_stoer
519 * \brief The Bulirsch-Stoer algorithm.
520 *
521 * The Bulirsch-Stoer is a controlled stepper that adjusts both step size
522 * and order of the method. The algorithm uses the modified midpoint and
523 * a polynomial extrapolation compute the solution.
524 *
525 * \tparam State The state type.
526 * \tparam Value The value type.
527 * \tparam Deriv The type representing the time derivative of the state.
528 * \tparam Time The time representing the independent variable - the time.
529 * \tparam Algebra The algebra type.
530 * \tparam Operations The operations type.
531 * \tparam Resizer The resizer policy type.
532 */
533
534 /**
535 * \fn bulirsch_stoer::bulirsch_stoer( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt )
536 * \brief Constructs the bulirsch_stoer class, including initialization of
537 * the error bounds.
538 *
539 * \param eps_abs Absolute tolerance level.
540 * \param eps_rel Relative tolerance level.
541 * \param factor_x Factor for the weight of the state.
542 * \param factor_dxdt Factor for the weight of the derivative.
543 */
544
545 /**
546 * \fn bulirsch_stoer::try_step( System system , StateInOut &x , time_type &t , time_type &dt )
547 * \brief Tries to perform one step.
548 *
549 * This method tries to do one step with step size dt. If the error estimate
550 * is to large, the step is rejected and the method returns fail and the
551 * step size dt is reduced. If the error estimate is acceptably small, the
552 * step is performed, success is returned and dt might be increased to make
553 * the steps as large as possible. This method also updates t if a step is
554 * performed. Also, the internal order of the stepper is adjusted if required.
555 *
556 * \param system The system function to solve, hence the r.h.s. of the ODE.
557 * It must fulfill the Simple System concept.
558 * \param x The state of the ODE which should be solved. Overwritten if
559 * the step is successful.
560 * \param t The value of the time. Updated if the step is successful.
561 * \param dt The step size. Updated.
562 * \return success if the step was accepted, fail otherwise.
563 */
564
565 /**
566 * \fn bulirsch_stoer::try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt )
567 * \brief Tries to perform one step.
568 *
569 * This method tries to do one step with step size dt. If the error estimate
570 * is to large, the step is rejected and the method returns fail and the
571 * step size dt is reduced. If the error estimate is acceptably small, the
572 * step is performed, success is returned and dt might be increased to make
573 * the steps as large as possible. This method also updates t if a step is
574 * performed. Also, the internal order of the stepper is adjusted if required.
575 *
576 * \param system The system function to solve, hence the r.h.s. of the ODE.
577 * It must fulfill the Simple System concept.
578 * \param x The state of the ODE which should be solved. Overwritten if
579 * the step is successful.
580 * \param dxdt The derivative of state.
581 * \param t The value of the time. Updated if the step is successful.
582 * \param dt The step size. Updated.
583 * \return success if the step was accepted, fail otherwise.
584 */
585
586 /**
587 * \fn bulirsch_stoer::try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt )
588 * \brief Tries to perform one step.
589 *
590 * \note This method is disabled if state_type=time_type to avoid ambiguity.
591 *
592 * This method tries to do one step with step size dt. If the error estimate
593 * is to large, the step is rejected and the method returns fail and the
594 * step size dt is reduced. If the error estimate is acceptably small, the
595 * step is performed, success is returned and dt might be increased to make
596 * the steps as large as possible. This method also updates t if a step is
597 * performed. Also, the internal order of the stepper is adjusted if required.
598 *
599 * \param system The system function to solve, hence the r.h.s. of the ODE.
600 * It must fulfill the Simple System concept.
601 * \param in The state of the ODE which should be solved.
602 * \param t The value of the time. Updated if the step is successful.
603 * \param out Used to store the result of the step.
604 * \param dt The step size. Updated.
605 * \return success if the step was accepted, fail otherwise.
606 */
607
608
609 /**
610 * \fn bulirsch_stoer::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
611 * \brief Tries to perform one step.
612 *
613 * This method tries to do one step with step size dt. If the error estimate
614 * is to large, the step is rejected and the method returns fail and the
615 * step size dt is reduced. If the error estimate is acceptably small, the
616 * step is performed, success is returned and dt might be increased to make
617 * the steps as large as possible. This method also updates t if a step is
618 * performed. Also, the internal order of the stepper is adjusted if required.
619 *
620 * \param system The system function to solve, hence the r.h.s. of the ODE.
621 * It must fulfill the Simple System concept.
622 * \param in The state of the ODE which should be solved.
623 * \param dxdt The derivative of state.
624 * \param t The value of the time. Updated if the step is successful.
625 * \param out Used to store the result of the step.
626 * \param dt The step size. Updated.
627 * \return success if the step was accepted, fail otherwise.
628 */
629
630
631 /**
632 * \fn bulirsch_stoer::adjust_size( const StateIn &x )
633 * \brief Adjust the size of all temporaries in the stepper manually.
634 * \param x A state from which the size of the temporaries to be resized is deduced.
635 */
636
637}
638}
639}
640
641#endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED
642

source code of boost/boost/numeric/odeint/stepper/bulirsch_stoer.hpp