| /* |
| [auto_generated] |
| libs/numeric/odeint/examples/black_hole.cpp |
| |
| [begin_description] |
| This example shows how the __float128 from gcc libquadmath can be used with odeint. |
| [end_description] |
| |
| Copyright 2012 Karsten Ahnert |
| Copyright 2012 Lee Hodgkinson |
| Copyright 2012 Mario Mulansky |
| |
| Distributed under the Boost Software License, Version 1.0. |
| (See accompanying file LICENSE_1_0.txt or |
| copy at http://www.boost.org/LICENSE_1_0.txt) |
| */ |
| |
| #include <cstdlib> |
| #include <cmath> |
| #include <iostream> |
| #include <iterator> |
| #include <utility> |
| #include <algorithm> |
| #include <cassert> |
| #include <vector> |
| #include <complex> |
| |
| extern "C" { |
| #include <quadmath.h> |
| } |
| |
| const __float128 zero =strtoflt128 ("0.0", NULL); |
| |
| namespace std { |
| |
| inline __float128 abs( __float128 x ) |
| { |
| return fabsq( x ); |
| } |
| |
| inline __float128 sqrt( __float128 x ) |
| { |
| return sqrtq( x ); |
| } |
| |
| inline __float128 pow( __float128 x , __float128 y ) |
| { |
| return powq( x , y ); |
| } |
| |
| inline __float128 abs( std::complex< __float128 > x ) |
| { |
| return sqrtq( x.real() * x.real() + x.imag() * x.imag() ); |
| } |
| |
| inline std::complex< __float128 > pow( std::complex< __float128> x , __float128 y ) |
| { |
| __float128 r = pow( abs(x) , y ); |
| __float128 phi = atanq( x.imag() / x.real() ); |
| return std::complex< __float128 >( r * cosq( y * phi ) , r * sinq( y * phi ) ); |
| } |
| } |
| |
| inline std::ostream& operator<< (std::ostream& os, const __float128& f) { |
| |
| char* y = new char[1000]; |
| quadmath_snprintf(y, 1000, "%.30Qg", f) ; |
| os.precision(30); |
| os<<y; |
| delete[] y; |
| return os; |
| } |
| |
| |
| #include <boost/array.hpp> |
| #include <boost/range/algorithm.hpp> |
| #include <boost/range/adaptor/filtered.hpp> |
| #include <boost/range/numeric.hpp> |
| #include <boost/numeric/odeint.hpp> |
| |
| |
| |
| using namespace boost::numeric::odeint; |
| using namespace std; |
| |
| typedef __float128 my_float; |
| typedef std::vector< std::complex < my_float > > state_type; |
| |
| struct radMod |
| { |
| my_float m_om; |
| my_float m_l; |
| |
| radMod( my_float om , my_float l ) |
| : m_om( om ) , m_l( l ) { } |
| |
| void operator()( const state_type &x , state_type &dxdt , my_float r ) const |
| { |
| |
| dxdt[0] = x[1]; |
| dxdt[1] = -(2*(r-1)/(r*(r-2)))*x[1]-((m_om*m_om*r*r/((r-2)*(r-2)))-(m_l*(m_l+1)/(r*(r-2))))*x[0]; |
| } |
| }; |
| |
| |
| |
| |
| |
| |
| |
| int main( int argc , char **argv ) |
| { |
| |
| |
| state_type x(2); |
| |
| my_float re0 = strtoflt128( "-0.00008944230755601224204687038354994353820468" , NULL ); |
| my_float im0 = strtoflt128( "0.00004472229441850588228136889483397204368247" , NULL ); |
| my_float re1 = strtoflt128( "-4.464175354293244250869336196695966076150E-6 " , NULL ); |
| my_float im1 = strtoflt128( "-8.950483248390306670770345406051469584488E-6" , NULL ); |
| |
| x[0] = complex< my_float >( re0 ,im0 ); |
| x[1] = complex< my_float >( re1 ,im1 ); |
| |
| const my_float dt =strtoflt128 ("-0.001", NULL); |
| const my_float start =strtoflt128 ("10000.0", NULL); |
| const my_float end =strtoflt128 ("9990.0", NULL); |
| const my_float omega =strtoflt128 ("2.0", NULL); |
| const my_float ell =strtoflt128 ("1.0", NULL); |
| |
| |
| |
| my_float abs_err = strtoflt128( "1.0E-15" , NULL ) , rel_err = strtoflt128( "1.0E-10" , NULL ); |
| my_float a_x = strtoflt128( "1.0" , NULL ) , a_dxdt = strtoflt128( "1.0" , NULL ); |
| |
| typedef runge_kutta_dopri5< state_type, my_float > dopri5_type; |
| typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type; |
| typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type; |
| |
| dense_output_dopri5_type dopri5( controlled_dopri5_type( default_error_checker< my_float >( abs_err , rel_err , a_x , a_dxdt ) ) ); |
| |
| std::for_each( make_adaptive_time_iterator_begin(dopri5 , radMod(omega , ell) , x , start , end , dt) , |
| make_adaptive_time_iterator_end(dopri5 , radMod(omega , ell) , x ) , |
| []( const std::pair< state_type&, my_float > &x ) { |
| std::cout << x.second << ", " << x.first[0].real() << "\n"; } |
| ); |
| |
| |
| |
| return 0; |
| } |