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 | |
49 | namespace boost { |
50 | namespace numeric { |
51 | namespace odeint { |
52 | |
53 | template< |
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 | > |
62 | class bulirsch_stoer { |
63 | |
64 | public: |
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 | |
362 | private: |
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 ( 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 | |