]>
Commit | Line | Data |
---|---|---|
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( 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 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 |