]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* |
2 | [auto_generated] | |
3 | boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp | |
4 | ||
5 | [begin_description] | |
6 | Implementaiton of the Burlish-Stoer method with dense output | |
7 | [end_description] | |
8 | ||
9 | Copyright 2011-2015 Mario Mulansky | |
10 | Copyright 2011-2013 Karsten Ahnert | |
11 | Copyright 2012 Christoph Koke | |
12 | ||
13 | Distributed under the Boost Software License, Version 1.0. | |
14 | (See accompanying file LICENSE_1_0.txt or | |
15 | copy at http://www.boost.org/LICENSE_1_0.txt) | |
16 | */ | |
17 | ||
18 | ||
19 | #ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED | |
20 | #define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_DENSE_OUT_HPP_INCLUDED | |
21 | ||
22 | ||
23 | #include <iostream> | |
24 | ||
25 | #include <algorithm> | |
26 | ||
27 | #include <boost/config.hpp> // for min/max guidelines | |
28 | ||
29 | #include <boost/numeric/odeint/util/bind.hpp> | |
30 | ||
31 | #include <boost/math/special_functions/binomial.hpp> | |
32 | ||
33 | #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp> | |
34 | #include <boost/numeric/odeint/stepper/modified_midpoint.hpp> | |
35 | #include <boost/numeric/odeint/stepper/controlled_step_result.hpp> | |
36 | #include <boost/numeric/odeint/algebra/range_algebra.hpp> | |
37 | #include <boost/numeric/odeint/algebra/default_operations.hpp> | |
38 | #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> | |
39 | #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> | |
40 | ||
41 | #include <boost/numeric/odeint/util/state_wrapper.hpp> | |
42 | #include <boost/numeric/odeint/util/is_resizeable.hpp> | |
43 | #include <boost/numeric/odeint/util/resizer.hpp> | |
44 | #include <boost/numeric/odeint/util/unit_helper.hpp> | |
45 | ||
46 | #include <boost/numeric/odeint/integrate/max_step_checker.hpp> | |
47 | ||
48 | #include <boost/type_traits.hpp> | |
49 | ||
50 | ||
51 | namespace boost { | |
52 | namespace numeric { | |
53 | namespace odeint { | |
54 | ||
55 | template< | |
56 | class State , | |
57 | class Value = double , | |
58 | class Deriv = State , | |
59 | class Time = Value , | |
60 | class Algebra = typename algebra_dispatcher< State >::algebra_type , | |
61 | class Operations = typename operations_dispatcher< State >::operations_type , | |
62 | class Resizer = initially_resizer | |
63 | > | |
64 | class bulirsch_stoer_dense_out { | |
65 | ||
66 | ||
67 | public: | |
68 | ||
69 | typedef State state_type; | |
70 | typedef Value value_type; | |
71 | typedef Deriv deriv_type; | |
72 | typedef Time time_type; | |
73 | typedef Algebra algebra_type; | |
74 | typedef Operations operations_type; | |
75 | typedef Resizer resizer_type; | |
76 | typedef dense_output_stepper_tag stepper_category; | |
77 | #ifndef DOXYGEN_SKIP | |
78 | typedef state_wrapper< state_type > wrapped_state_type; | |
79 | typedef state_wrapper< deriv_type > wrapped_deriv_type; | |
80 | ||
81 | typedef bulirsch_stoer_dense_out< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type; | |
82 | ||
83 | typedef typename inverse_time< time_type >::type inv_time_type; | |
84 | ||
85 | typedef std::vector< value_type > value_vector; | |
86 | typedef std::vector< time_type > time_vector; | |
87 | typedef std::vector< inv_time_type > inv_time_vector; //should be 1/time_type for boost.units | |
88 | typedef std::vector< value_vector > value_matrix; | |
89 | typedef std::vector< size_t > int_vector; | |
90 | typedef std::vector< wrapped_state_type > state_vector_type; | |
91 | typedef std::vector< wrapped_deriv_type > deriv_vector_type; | |
92 | typedef std::vector< deriv_vector_type > deriv_table_type; | |
93 | #endif //DOXYGEN_SKIP | |
94 | ||
95 | const static size_t m_k_max = 8; | |
96 | ||
97 | ||
98 | ||
99 | bulirsch_stoer_dense_out( | |
100 | value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 , | |
101 | value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 , | |
102 | time_type max_dt = static_cast<time_type>(0) , | |
103 | bool control_interpolation = false ) | |
104 | : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , | |
105 | m_max_dt(max_dt) , | |
106 | m_control_interpolation( control_interpolation) , | |
107 | m_last_step_rejected( false ) , m_first( true ) , | |
108 | m_current_state_x1( true ) , | |
109 | m_error( m_k_max ) , | |
110 | m_interval_sequence( m_k_max+1 ) , | |
111 | m_coeff( m_k_max+1 ) , | |
112 | m_cost( m_k_max+1 ) , | |
b32b8144 | 113 | m_facmin_table( m_k_max+1 ) , |
7c673cae FG |
114 | m_table( m_k_max ) , |
115 | m_mp_states( m_k_max+1 ) , | |
116 | m_derivs( m_k_max+1 ) , | |
117 | m_diffs( 2*m_k_max+2 ) , | |
118 | STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 ) | |
119 | { | |
120 | BOOST_USING_STD_MIN(); | |
121 | BOOST_USING_STD_MAX(); | |
122 | ||
123 | for( unsigned short i = 0; i < m_k_max+1; i++ ) | |
124 | { | |
125 | /* only this specific sequence allows for dense output */ | |
126 | m_interval_sequence[i] = 2 + 4*i; // 2 6 10 14 ... | |
127 | m_derivs[i].resize( m_interval_sequence[i] ); | |
128 | if( i == 0 ) | |
b32b8144 | 129 | { |
7c673cae | 130 | m_cost[i] = m_interval_sequence[i]; |
b32b8144 FG |
131 | } else |
132 | { | |
7c673cae | 133 | m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; |
b32b8144 FG |
134 | } |
135 | m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) ); | |
7c673cae FG |
136 | m_coeff[i].resize(i); |
137 | for( size_t k = 0 ; k < i ; ++k ) | |
138 | { | |
139 | const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] ); | |
140 | m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation | |
141 | } | |
142 | // crude estimate of optimal order | |
143 | ||
144 | m_current_k_opt = 4; | |
145 | /* no calculation because log10 might not exist for value_type! | |
146 | const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >( 1.0E-12 ) ) ) * 0.6 + 0.5 ); | |
147 | m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 1 , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>( m_k_max-1 ) , static_cast<int>( logfact ) )); | |
148 | */ | |
149 | } | |
150 | int num = 1; | |
151 | for( int i = 2*(m_k_max)+1 ; i >=0 ; i-- ) | |
152 | { | |
153 | m_diffs[i].resize( num ); | |
154 | num += (i+1)%2; | |
155 | } | |
156 | } | |
157 | ||
158 | template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut > | |
159 | controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt ) | |
160 | { | |
161 | if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) ) | |
162 | { | |
163 | // given step size is bigger then max_dt | |
164 | // set limit and return fail | |
165 | dt = m_max_dt; | |
166 | return fail; | |
167 | } | |
168 | ||
169 | BOOST_USING_STD_MIN(); | |
170 | BOOST_USING_STD_MAX(); | |
171 | using std::pow; | |
172 | ||
173 | static const value_type val1( 1.0 ); | |
174 | ||
175 | bool reject( true ); | |
176 | ||
177 | time_vector h_opt( m_k_max+1 ); | |
178 | inv_time_vector work( m_k_max+1 ); | |
179 | ||
180 | m_k_final = 0; | |
181 | time_type new_h = dt; | |
182 | ||
183 | //std::cout << "t=" << t <<", dt=" << dt << ", k_opt=" << m_current_k_opt << ", first: " << m_first << std::endl; | |
184 | ||
185 | for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ ) | |
186 | { | |
187 | m_midpoint.set_steps( m_interval_sequence[k] ); | |
188 | if( k == 0 ) | |
189 | { | |
190 | m_midpoint.do_step( system , in , dxdt , t , out , dt , m_mp_states[k].m_v , m_derivs[k]); | |
191 | } | |
192 | else | |
193 | { | |
194 | m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt , m_mp_states[k].m_v , m_derivs[k] ); | |
195 | extrapolate( k , m_table , m_coeff , out ); | |
196 | // get error estimate | |
197 | m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v , | |
198 | typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) ); | |
199 | const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt ); | |
200 | h_opt[k] = calc_h_opt( dt , error , k ); | |
201 | work[k] = static_cast<value_type>( m_cost[k] ) / h_opt[k]; | |
202 | ||
203 | m_k_final = k; | |
204 | ||
205 | if( (k == m_current_k_opt-1) || m_first ) | |
206 | { // convergence before k_opt ? | |
207 | if( error < 1.0 ) | |
208 | { | |
209 | //convergence | |
210 | reject = false; | |
211 | if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) ) | |
212 | { | |
213 | // leave order as is (except we were in first round) | |
214 | 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 ) ); | |
215 | new_h = h_opt[k] * static_cast<value_type>( m_cost[k+1] ) / static_cast<value_type>( m_cost[k] ); | |
216 | } else { | |
217 | 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) ) ); | |
218 | new_h = h_opt[k]; | |
219 | } | |
220 | break; | |
221 | } | |
222 | else if( should_reject( error , k ) && !m_first ) | |
223 | { | |
224 | reject = true; | |
225 | new_h = h_opt[k]; | |
226 | break; | |
227 | } | |
228 | } | |
229 | if( k == m_current_k_opt ) | |
230 | { // convergence at k_opt ? | |
231 | if( error < 1.0 ) | |
232 | { | |
233 | //convergence | |
234 | reject = false; | |
235 | if( (work[k-1] < KFAC2*work[k]) ) | |
236 | { | |
237 | m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 ); | |
238 | new_h = h_opt[m_current_k_opt]; | |
239 | } | |
240 | else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected ) | |
241 | { | |
242 | m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(m_current_k_opt)+1 ); | |
243 | new_h = h_opt[k]*static_cast<value_type>( m_cost[m_current_k_opt] ) / static_cast<value_type>( m_cost[k] ); | |
244 | } else | |
245 | new_h = h_opt[m_current_k_opt]; | |
246 | break; | |
247 | } | |
248 | else if( should_reject( error , k ) ) | |
249 | { | |
250 | reject = true; | |
251 | new_h = h_opt[m_current_k_opt]; | |
252 | break; | |
253 | } | |
254 | } | |
255 | if( k == m_current_k_opt+1 ) | |
256 | { // convergence at k_opt+1 ? | |
257 | if( error < 1.0 ) | |
258 | { //convergence | |
259 | reject = false; | |
260 | if( work[k-2] < KFAC2*work[k-1] ) | |
261 | m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast<int>(m_current_k_opt)-1 ); | |
262 | if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected ) | |
263 | m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<int>(m_k_max)-1 , static_cast<int>(k) ); | |
264 | new_h = h_opt[m_current_k_opt]; | |
265 | } else | |
266 | { | |
267 | reject = true; | |
268 | new_h = h_opt[m_current_k_opt]; | |
269 | } | |
270 | break; | |
271 | } | |
272 | } | |
273 | } | |
274 | ||
275 | if( !reject ) | |
276 | { | |
277 | ||
278 | //calculate dxdt for next step and dense output | |
279 | typename odeint::unwrap_reference< System >::type &sys = system; | |
280 | sys( out , dxdt_new , t+dt ); | |
281 | ||
282 | //prepare dense output | |
283 | value_type error = prepare_dense_output( m_k_final , in , dxdt , out , dxdt_new , dt ); | |
284 | ||
285 | if( error > static_cast<value_type>(10) ) // we are not as accurate for interpolation as for the steps | |
286 | { | |
287 | reject = true; | |
288 | new_h = dt * pow BOOST_PREVENT_MACRO_SUBSTITUTION( error , static_cast<value_type>(-1)/(2*m_k_final+2) ); | |
289 | } else { | |
290 | t += dt; | |
291 | } | |
292 | } | |
293 | //set next stepsize | |
294 | if( !m_last_step_rejected || (new_h < dt) ) | |
295 | { | |
296 | // limit step size | |
297 | if( m_max_dt != static_cast<time_type>(0) ) | |
298 | { | |
299 | new_h = detail::min_abs(m_max_dt, new_h); | |
300 | } | |
301 | dt = new_h; | |
302 | } | |
303 | ||
304 | m_last_step_rejected = reject; | |
305 | if( reject ) | |
306 | return fail; | |
307 | else | |
308 | return success; | |
309 | } | |
310 | ||
311 | template< class StateType > | |
312 | void initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 ) | |
313 | { | |
314 | m_resizer.adjust_size( x0 , detail::bind( &controlled_error_bs_type::template resize_impl< StateType > , detail::ref( *this ) , detail::_1 ) ); | |
315 | boost::numeric::odeint::copy( x0 , get_current_state() ); | |
316 | m_t = t0; | |
317 | m_dt = dt0; | |
318 | reset(); | |
319 | } | |
320 | ||
321 | ||
322 | /* ======================================================= | |
323 | * the actual step method that should be called from outside (maybe make try_step private?) | |
324 | */ | |
325 | template< class System > | |
326 | std::pair< time_type , time_type > do_step( System system ) | |
327 | { | |
328 | if( m_first ) | |
329 | { | |
330 | typename odeint::unwrap_reference< System >::type &sys = system; | |
331 | sys( get_current_state() , get_current_deriv() , m_t ); | |
332 | } | |
333 | ||
334 | failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails | |
335 | controlled_step_result res = fail; | |
336 | m_t_last = m_t; | |
337 | while( res == fail ) | |
338 | { | |
339 | res = try_step( system , get_current_state() , get_current_deriv() , m_t , get_old_state() , get_old_deriv() , m_dt ); | |
340 | m_first = false; | |
341 | fail_checker(); // check for overflow of failed steps | |
342 | } | |
343 | toggle_current_state(); | |
344 | return std::make_pair( m_t_last , m_t ); | |
345 | } | |
346 | ||
347 | /* performs the interpolation from a calculated step */ | |
348 | template< class StateOut > | |
349 | void calc_state( time_type t , StateOut &x ) const | |
350 | { | |
351 | do_interpolation( t , x ); | |
352 | } | |
353 | ||
354 | const state_type& current_state( void ) const | |
355 | { | |
356 | return get_current_state(); | |
357 | } | |
358 | ||
359 | time_type current_time( void ) const | |
360 | { | |
361 | return m_t; | |
362 | } | |
363 | ||
364 | const state_type& previous_state( void ) const | |
365 | { | |
366 | return get_old_state(); | |
367 | } | |
368 | ||
369 | time_type previous_time( void ) const | |
370 | { | |
371 | return m_t_last; | |
372 | } | |
373 | ||
374 | time_type current_time_step( void ) const | |
375 | { | |
376 | return m_dt; | |
377 | } | |
378 | ||
379 | /** \brief Resets the internal state of the stepper. */ | |
380 | void reset() | |
381 | { | |
382 | m_first = true; | |
383 | m_last_step_rejected = false; | |
384 | } | |
385 | ||
386 | template< class StateIn > | |
387 | void adjust_size( const StateIn &x ) | |
388 | { | |
389 | resize_impl( x ); | |
390 | m_midpoint.adjust_size( x ); | |
391 | } | |
392 | ||
393 | ||
92f5a8d4 TL |
394 | protected: |
395 | ||
396 | time_type m_max_dt; | |
397 | ||
398 | ||
7c673cae FG |
399 | private: |
400 | ||
401 | template< class StateInOut , class StateVector > | |
402 | void extrapolate( size_t k , StateVector &table , const value_matrix &coeff , StateInOut &xest , size_t order_start_index = 0 ) | |
403 | //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf | |
404 | { | |
405 | static const value_type val1( 1.0 ); | |
406 | for( int j=k-1 ; j>0 ; --j ) | |
407 | { | |
408 | m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v , | |
409 | typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index] , | |
410 | -coeff[k + order_start_index][j + order_start_index] ) ); | |
411 | } | |
412 | m_algebra.for_each3( xest , table[0].m_v , xest , | |
413 | typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][0 + order_start_index] , | |
414 | -coeff[k + order_start_index][0 + order_start_index]) ); | |
415 | } | |
416 | ||
417 | ||
418 | template< class StateVector > | |
419 | void extrapolate_dense_out( size_t k , StateVector &table , const value_matrix &coeff , size_t order_start_index = 0 ) | |
420 | //polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf | |
421 | { | |
422 | // result is written into table[0] | |
423 | static const value_type val1( 1.0 ); | |
424 | for( int j=k ; j>1 ; --j ) | |
425 | { | |
426 | m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v , | |
427 | typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][j + order_start_index - 1] , | |
428 | -coeff[k + order_start_index][j + order_start_index - 1] ) ); | |
429 | } | |
430 | m_algebra.for_each3( table[0].m_v , table[1].m_v , table[0].m_v , | |
431 | typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k + order_start_index][order_start_index] , | |
432 | -coeff[k + order_start_index][order_start_index]) ); | |
433 | } | |
434 | ||
435 | time_type calc_h_opt( time_type h , value_type error , size_t k ) const | |
436 | { | |
437 | BOOST_USING_STD_MIN(); | |
438 | BOOST_USING_STD_MAX(); | |
439 | using std::pow; | |
440 | ||
441 | value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]); | |
b32b8144 | 442 | value_type facmin = m_facmin_table[k]; |
7c673cae FG |
443 | value_type fac; |
444 | if (error == 0.0) | |
445 | fac = static_cast<value_type>(1)/facmin; | |
446 | else | |
447 | { | |
448 | fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo ); | |
449 | fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( facmin/STEPFAC4 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(static_cast<value_type>(1)/facmin) , fac ) ); | |
450 | } | |
451 | return h*fac; | |
452 | } | |
453 | ||
454 | bool in_convergence_window( size_t k ) const | |
455 | { | |
456 | if( (k == m_current_k_opt-1) && !m_last_step_rejected ) | |
457 | return true; // decrease order only if last step was not rejected | |
458 | return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) ); | |
459 | } | |
460 | ||
461 | bool should_reject( value_type error , size_t k ) const | |
462 | { | |
463 | if( k == m_current_k_opt-1 ) | |
464 | { | |
465 | const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] / | |
466 | (m_interval_sequence[0]*m_interval_sequence[0]); | |
467 | //step will fail, criterion 17.3.17 in NR | |
468 | return ( error > d*d ); | |
469 | } | |
470 | else if( k == m_current_k_opt ) | |
471 | { | |
472 | const value_type d = m_interval_sequence[m_current_k_opt+1] / m_interval_sequence[0]; | |
473 | return ( error > d*d ); | |
474 | } else | |
475 | return error > 1.0; | |
476 | } | |
477 | ||
478 | template< class StateIn1 , class DerivIn1 , class StateIn2 , class DerivIn2 > | |
479 | value_type prepare_dense_output( int k , const StateIn1 &x_start , const DerivIn1 &dxdt_start , | |
480 | const StateIn2 & /* x_end */ , const DerivIn2 & /*dxdt_end */ , time_type dt ) | |
481 | /* k is the order to which the result was approximated */ | |
482 | { | |
483 | ||
484 | /* compute the coefficients of the interpolation polynomial | |
485 | * we parametrize the interval t .. t+dt by theta = -1 .. 1 | |
486 | * we use 2k+3 values at the interval center theta=0 to obtain the interpolation coefficients | |
487 | * the values are x(t+dt/2) and the derivatives dx/dt , ... d^(2k+2) x / dt^(2k+2) at the midpoints | |
488 | * the derivatives are approximated via finite differences | |
489 | * all values are obtained from interpolation of the results from the increasing orders of the midpoint calls | |
490 | */ | |
491 | ||
492 | // calculate finite difference approximations to derivatives at the midpoint | |
493 | for( int j = 0 ; j<=k ; j++ ) | |
494 | { | |
495 | /* not working with boost units... */ | |
496 | const value_type d = m_interval_sequence[j] / ( static_cast<value_type>(2) * dt ); | |
497 | value_type f = 1.0; //factor 1/2 here because our interpolation interval has length 2 !!! | |
498 | for( int kappa = 0 ; kappa <= 2*j+1 ; ++kappa ) | |
499 | { | |
500 | calculate_finite_difference( j , kappa , f , dxdt_start ); | |
501 | f *= d; | |
502 | } | |
503 | ||
504 | if( j > 0 ) | |
505 | extrapolate_dense_out( j , m_mp_states , m_coeff ); | |
506 | } | |
507 | ||
508 | time_type d = dt/2; | |
509 | ||
510 | // extrapolate finite differences | |
511 | for( int kappa = 0 ; kappa<=2*k+1 ; kappa++ ) | |
512 | { | |
513 | for( int j=1 ; j<=(k-kappa/2) ; ++j ) | |
514 | extrapolate_dense_out( j , m_diffs[kappa] , m_coeff , kappa/2 ); | |
515 | ||
516 | // extrapolation results are now stored in m_diffs[kappa][0] | |
517 | ||
518 | // divide kappa-th derivative by kappa because we need these terms for dense output interpolation | |
519 | m_algebra.for_each1( m_diffs[kappa][0].m_v , typename operations_type::template scale< time_type >( static_cast<time_type>(d) ) ); | |
520 | ||
521 | d *= dt/(2*(kappa+2)); | |
522 | } | |
523 | ||
524 | // dense output coefficients a_0 is stored in m_mp_states[0], a_i for i = 1...2k are stored in m_diffs[i-1][0] | |
525 | ||
526 | // the error is just the highest order coefficient of the interpolation polynomial | |
527 | // this is because we use only the midpoint theta=0 as support for the interpolation (remember that theta = -1 .. 1) | |
528 | ||
529 | value_type error = 0.0; | |
530 | if( m_control_interpolation ) | |
531 | { | |
532 | boost::numeric::odeint::copy( m_diffs[2*k+1][0].m_v , m_err.m_v ); | |
533 | error = m_error_checker.error( m_algebra , x_start , dxdt_start , m_err.m_v , dt ); | |
534 | } | |
535 | ||
536 | return error; | |
537 | } | |
538 | ||
539 | template< class DerivIn > | |
540 | void calculate_finite_difference( size_t j , size_t kappa , value_type fac , const DerivIn &dxdt ) | |
541 | { | |
542 | const int m = m_interval_sequence[j]/2-1; | |
543 | if( kappa == 0) // no calculation required for 0th derivative of f | |
544 | { | |
545 | m_algebra.for_each2( m_diffs[0][j].m_v , m_derivs[j][m].m_v , | |
546 | typename operations_type::template scale_sum1< value_type >( fac ) ); | |
547 | } | |
548 | else | |
549 | { | |
550 | // calculate the index of m_diffs for this kappa-j-combination | |
551 | const int j_diffs = j - kappa/2; | |
552 | ||
553 | m_algebra.for_each2( m_diffs[kappa][j_diffs].m_v , m_derivs[j][m+kappa].m_v , | |
554 | typename operations_type::template scale_sum1< value_type >( fac ) ); | |
555 | value_type sign = -1.0; | |
556 | int c = 1; | |
557 | //computes the j-th order finite difference for the kappa-th derivative of f at t+dt/2 using function evaluations stored in m_derivs | |
558 | for( int i = m+static_cast<int>(kappa)-2 ; i >= m-static_cast<int>(kappa) ; i -= 2 ) | |
559 | { | |
560 | if( i >= 0 ) | |
561 | { | |
562 | m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , m_derivs[j][i].m_v , | |
563 | typename operations_type::template scale_sum2< value_type , value_type >( 1.0 , | |
564 | sign * fac * boost::math::binomial_coefficient< value_type >( kappa , c ) ) ); | |
565 | } | |
566 | else | |
567 | { | |
568 | m_algebra.for_each3( m_diffs[kappa][j_diffs].m_v , m_diffs[kappa][j_diffs].m_v , dxdt , | |
569 | typename operations_type::template scale_sum2< value_type , value_type >( 1.0 , sign * fac ) ); | |
570 | } | |
571 | sign *= -1; | |
572 | ++c; | |
573 | } | |
574 | } | |
575 | } | |
576 | ||
577 | template< class StateOut > | |
578 | void do_interpolation( time_type t , StateOut &out ) const | |
579 | { | |
580 | // interpolation polynomial is defined for theta = -1 ... 1 | |
581 | // m_k_final is the number of order-iterations done for the last step - it governs the order of the interpolation polynomial | |
582 | const value_type theta = 2 * get_unit_value( (t - m_t_last) / (m_t - m_t_last) ) - 1; | |
583 | // we use only values at interval center, that is theta=0, for interpolation | |
584 | // our interpolation polynomial is thus of order 2k+2, hence we have 2k+3 terms | |
585 | ||
586 | boost::numeric::odeint::copy( m_mp_states[0].m_v , out ); | |
587 | // add remaining terms: x += a_1 theta + a2 theta^2 + ... + a_{2k} theta^{2k} | |
588 | value_type theta_pow( theta ); | |
589 | for( size_t i=0 ; i<=2*m_k_final+1 ; ++i ) | |
590 | { | |
591 | m_algebra.for_each3( out , out , m_diffs[i][0].m_v , | |
592 | typename operations_type::template scale_sum2< value_type >( static_cast<value_type>(1) , theta_pow ) ); | |
593 | theta_pow *= theta; | |
594 | } | |
595 | } | |
596 | ||
597 | /* Resizer methods */ | |
598 | template< class StateIn > | |
599 | bool resize_impl( const StateIn &x ) | |
600 | { | |
601 | bool resized( false ); | |
602 | ||
603 | resized |= adjust_size_by_resizeability( m_x1 , x , typename is_resizeable<state_type>::type() ); | |
604 | resized |= adjust_size_by_resizeability( m_x2 , x , typename is_resizeable<state_type>::type() ); | |
605 | resized |= adjust_size_by_resizeability( m_dxdt1 , x , typename is_resizeable<state_type>::type() ); | |
606 | resized |= adjust_size_by_resizeability( m_dxdt2 , x , typename is_resizeable<state_type>::type() ); | |
607 | resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable<state_type>::type() ); | |
608 | ||
609 | for( size_t i = 0 ; i < m_k_max ; ++i ) | |
610 | resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable<state_type>::type() ); | |
611 | for( size_t i = 0 ; i < m_k_max+1 ; ++i ) | |
612 | resized |= adjust_size_by_resizeability( m_mp_states[i] , x , typename is_resizeable<state_type>::type() ); | |
613 | for( size_t i = 0 ; i < m_k_max+1 ; ++i ) | |
614 | for( size_t j = 0 ; j < m_derivs[i].size() ; ++j ) | |
615 | resized |= adjust_size_by_resizeability( m_derivs[i][j] , x , typename is_resizeable<deriv_type>::type() ); | |
616 | for( size_t i = 0 ; i < 2*m_k_max+2 ; ++i ) | |
617 | for( size_t j = 0 ; j < m_diffs[i].size() ; ++j ) | |
618 | resized |= adjust_size_by_resizeability( m_diffs[i][j] , x , typename is_resizeable<deriv_type>::type() ); | |
619 | ||
620 | return resized; | |
621 | } | |
622 | ||
623 | ||
624 | state_type& get_current_state( void ) | |
625 | { | |
626 | return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; | |
627 | } | |
628 | ||
629 | const state_type& get_current_state( void ) const | |
630 | { | |
631 | return m_current_state_x1 ? m_x1.m_v : m_x2.m_v ; | |
632 | } | |
633 | ||
634 | state_type& get_old_state( void ) | |
635 | { | |
636 | return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; | |
637 | } | |
638 | ||
639 | const state_type& get_old_state( void ) const | |
640 | { | |
641 | return m_current_state_x1 ? m_x2.m_v : m_x1.m_v ; | |
642 | } | |
643 | ||
644 | deriv_type& get_current_deriv( void ) | |
645 | { | |
646 | return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ; | |
647 | } | |
648 | ||
649 | const deriv_type& get_current_deriv( void ) const | |
650 | { | |
651 | return m_current_state_x1 ? m_dxdt1.m_v : m_dxdt2.m_v ; | |
652 | } | |
653 | ||
654 | deriv_type& get_old_deriv( void ) | |
655 | { | |
656 | return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ; | |
657 | } | |
658 | ||
659 | const deriv_type& get_old_deriv( void ) const | |
660 | { | |
661 | return m_current_state_x1 ? m_dxdt2.m_v : m_dxdt1.m_v ; | |
662 | } | |
663 | ||
664 | ||
665 | void toggle_current_state( void ) | |
666 | { | |
667 | m_current_state_x1 = ! m_current_state_x1; | |
668 | } | |
669 | ||
670 | ||
671 | ||
672 | default_error_checker< value_type, algebra_type , operations_type > m_error_checker; | |
673 | modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint; | |
674 | ||
7c673cae FG |
675 | bool m_control_interpolation; |
676 | ||
677 | bool m_last_step_rejected; | |
678 | bool m_first; | |
679 | ||
680 | time_type m_t; | |
681 | time_type m_dt; | |
682 | time_type m_dt_last; | |
683 | time_type m_t_last; | |
684 | ||
685 | size_t m_current_k_opt; | |
686 | size_t m_k_final; | |
687 | ||
688 | algebra_type m_algebra; | |
689 | ||
690 | resizer_type m_resizer; | |
691 | ||
692 | wrapped_state_type m_x1 , m_x2; | |
693 | wrapped_deriv_type m_dxdt1 , m_dxdt2; | |
694 | wrapped_state_type m_err; | |
695 | bool m_current_state_x1; | |
696 | ||
697 | ||
698 | ||
699 | value_vector m_error; // errors of repeated midpoint steps and extrapolations | |
700 | int_vector m_interval_sequence; // stores the successive interval counts | |
701 | value_matrix m_coeff; | |
702 | int_vector m_cost; // costs for interval count | |
b32b8144 | 703 | value_vector m_facmin_table; // for precomputed facmin to save pow calls |
7c673cae FG |
704 | |
705 | state_vector_type m_table; // sequence of states for extrapolation | |
706 | ||
707 | //for dense output: | |
708 | state_vector_type m_mp_states; // sequence of approximations of x at distance center | |
709 | deriv_table_type m_derivs; // table of function values | |
710 | deriv_table_type m_diffs; // table of function values | |
711 | ||
712 | //wrapped_state_type m_a1 , m_a2 , m_a3 , m_a4; | |
713 | ||
714 | value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2; | |
715 | }; | |
716 | ||
717 | ||
718 | ||
719 | /********** DOXYGEN **********/ | |
720 | ||
721 | /** | |
722 | * \class bulirsch_stoer_dense_out | |
723 | * \brief The Bulirsch-Stoer algorithm. | |
724 | * | |
725 | * The Bulirsch-Stoer is a controlled stepper that adjusts both step size | |
726 | * and order of the method. The algorithm uses the modified midpoint and | |
727 | * a polynomial extrapolation compute the solution. This class also provides | |
728 | * dense output facility. | |
729 | * | |
730 | * \tparam State The state type. | |
731 | * \tparam Value The value type. | |
732 | * \tparam Deriv The type representing the time derivative of the state. | |
733 | * \tparam Time The time representing the independent variable - the time. | |
734 | * \tparam Algebra The algebra type. | |
735 | * \tparam Operations The operations type. | |
736 | * \tparam Resizer The resizer policy type. | |
737 | */ | |
738 | ||
739 | /** | |
740 | * \fn bulirsch_stoer_dense_out::bulirsch_stoer_dense_out( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt , bool control_interpolation ) | |
741 | * \brief Constructs the bulirsch_stoer class, including initialization of | |
742 | * the error bounds. | |
743 | * | |
744 | * \param eps_abs Absolute tolerance level. | |
745 | * \param eps_rel Relative tolerance level. | |
746 | * \param factor_x Factor for the weight of the state. | |
747 | * \param factor_dxdt Factor for the weight of the derivative. | |
748 | * \param control_interpolation Set true to additionally control the error of | |
749 | * the interpolation. | |
750 | */ | |
751 | ||
752 | /** | |
753 | * \fn bulirsch_stoer_dense_out::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt ) | |
754 | * \brief Tries to perform one step. | |
755 | * | |
756 | * This method tries to do one step with step size dt. If the error estimate | |
757 | * is to large, the step is rejected and the method returns fail and the | |
758 | * step size dt is reduced. If the error estimate is acceptably small, the | |
759 | * step is performed, success is returned and dt might be increased to make | |
760 | * the steps as large as possible. This method also updates t if a step is | |
761 | * performed. Also, the internal order of the stepper is adjusted if required. | |
762 | * | |
763 | * \param system The system function to solve, hence the r.h.s. of the ODE. | |
764 | * It must fulfill the Simple System concept. | |
765 | * \param in The state of the ODE which should be solved. | |
766 | * \param dxdt The derivative of state. | |
767 | * \param t The value of the time. Updated if the step is successful. | |
768 | * \param out Used to store the result of the step. | |
769 | * \param dt The step size. Updated. | |
770 | * \return success if the step was accepted, fail otherwise. | |
771 | */ | |
772 | ||
773 | /** | |
774 | * \fn bulirsch_stoer_dense_out::initialize( const StateType &x0 , const time_type &t0 , const time_type &dt0 ) | |
775 | * \brief Initializes the dense output stepper. | |
776 | * | |
777 | * \param x0 The initial state. | |
778 | * \param t0 The initial time. | |
779 | * \param dt0 The initial time step. | |
780 | */ | |
781 | ||
782 | /** | |
783 | * \fn bulirsch_stoer_dense_out::do_step( System system ) | |
784 | * \brief Does one time step. This is the main method that should be used to | |
785 | * integrate an ODE with this stepper. | |
786 | * \note initialize has to be called before using this method to set the | |
787 | * initial conditions x,t and the stepsize. | |
788 | * \param system The system function to solve, hence the r.h.s. of the | |
789 | * ordinary differential equation. It must fulfill the Simple System concept. | |
790 | * \return Pair with start and end time of the integration step. | |
791 | */ | |
792 | ||
793 | /** | |
794 | * \fn bulirsch_stoer_dense_out::calc_state( time_type t , StateOut &x ) const | |
795 | * \brief Calculates the solution at an intermediate point within the last step | |
796 | * \param t The time at which the solution should be calculated, has to be | |
797 | * in the current time interval. | |
798 | * \param x The output variable where the result is written into. | |
799 | */ | |
800 | ||
801 | /** | |
802 | * \fn bulirsch_stoer_dense_out::current_state( void ) const | |
803 | * \brief Returns the current state of the solution. | |
804 | * \return The current state of the solution x(t). | |
805 | */ | |
806 | ||
807 | /** | |
808 | * \fn bulirsch_stoer_dense_out::current_time( void ) const | |
809 | * \brief Returns the current time of the solution. | |
810 | * \return The current time of the solution t. | |
811 | */ | |
812 | ||
813 | /** | |
814 | * \fn bulirsch_stoer_dense_out::previous_state( void ) const | |
815 | * \brief Returns the last state of the solution. | |
816 | * \return The last state of the solution x(t-dt). | |
817 | */ | |
818 | ||
819 | /** | |
820 | * \fn bulirsch_stoer_dense_out::previous_time( void ) const | |
821 | * \brief Returns the last time of the solution. | |
822 | * \return The last time of the solution t-dt. | |
823 | */ | |
824 | ||
825 | /** | |
826 | * \fn bulirsch_stoer_dense_out::current_time_step( void ) const | |
827 | * \brief Returns the current step size. | |
828 | * \return The current step size. | |
829 | */ | |
830 | ||
831 | /** | |
832 | * \fn bulirsch_stoer_dense_out::adjust_size( const StateIn &x ) | |
833 | * \brief Adjust the size of all temporaries in the stepper manually. | |
834 | * \param x A state from which the size of the temporaries to be resized is deduced. | |
835 | */ | |
836 | ||
837 | } | |
838 | } | |
839 | } | |
840 | ||
841 | #endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED |