gistfile1.cpp

#include <iostream>
#include <fstream>
#include <complex>
#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <mpfr.h>
#include <gmpxx.h>
#include "mpreal.h"

#include <boost/numeric/odeint.hpp>





using namespace std;
using namespace boost::numeric::odeint;
using mpfr::mpreal;


typedef mpreal my_float; //will set this to arb precision later.
typedef boost::numeric::ublas::vector<my_float > vector_type;
typedef boost::numeric::ublas::matrix<my_float > matrix_type;

//globals, will be used in order to keep precision high.
const my_float one = "1.0";
const my_float two = "2.0";

//[ radMod_system_function
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 vector_type &x , vector_type &dxdr , my_float r ) const
    {
        //x[0]=x0R,x[1]=x0I, x[2]=x1R, x[3]=x1I (leads to two sets of decoupled eqns)
        dxdr[0] = x[2];
        dxdr[2] = -(two*(r-one)/(r*(r-two)))*x[2]-((m_om*m_om*r*r/((r-two)*(r-two)))-(m_l*(m_l+one)/(r*(r-two))))*x[0];

        dxdr[1] = x[3];
        dxdr[3] = -(two*(r-one)/(r*(r-two)))*x[3]-((m_om*m_om*r*r/((r-two)*(r-two)))-(m_l*(m_l+one)/(r*(r-two))))*x[1];
    }
};
//]


struct radMod_jacobi
{

    my_float m_om;
    my_float m_l;

    radMod_jacobi( my_float om , my_float l) : m_om( om ) , m_l( l ) { }

    void operator()( const vector_type & x , matrix_type &J , const my_float & r , vector_type &dfdt )
    {
        
        J( 0 , 0 ) = 0;
        J( 0 , 1 ) = 0;
        J( 0 , 2 ) = 1;
        J( 0 , 3 ) = 0;

        J( 1 , 0 ) = 0;
        J( 1 , 1 ) = 0;
        J( 1 , 2 ) = 0;
        J( 1 , 3 ) = 1;

        J( 2 , 0 ) = -((m_om*m_om*r*r/((r-two)*(r-two)))-(m_l*(m_l+one)/(r*(r-two))));
        J( 2 , 1 ) = 0;
        J( 2 , 2 ) = -(two*(r-one)/(r*(r-two)));
        J( 2 , 3 ) = 0;

        J( 3 , 0 ) = 0;
        J( 3 , 1 ) = -((m_om*m_om*r*r/((r-two)*(r-two)))-(m_l*(m_l+one)/(r*(r-two))));
        J( 3 , 2 ) = 0;
        J( 3 , 3 ) = -(two*(r-one)/(r*(r-two)));
        



        dfdt[0] = 0.0;
        dfdt[1] = 0.0;
        dfdt[2] = two*(two-two*r+r*r)/(r*r*(r-two)*(r-two))*x[2]+((two*two*r*m_om*m_om/((r-two)*(r-two)*(r-two)))-(two*m_l*(m_l+one)*(r-one)/(r*r*(r-two)*(r-two))))*x[0];
        dfdt[3] = two*(two-two*r+r*r)/(r*r*(r-two)*(r-two))*x[3]+((two*two*r*m_om*m_om/((r-two)*(r-two)*(r-two)))-(two*m_l*(m_l+one)*(r-one)/(r*r*(r-two)*(r-two))))*x[1];
    }
};

struct streaming_observer
{
    std::ofstream& m_yout;
    std::ofstream& m_ypout;

    streaming_observer( std::ofstream &yout, std::ofstream &ypout ) : m_yout( yout ),m_ypout( ypout ) { }

    template < class State >
    void operator()( const std::pair< State& , my_float > &x ) const
    {
        // if(x.second > 10 && x.second <100){ //only output small range.

        m_yout.precision(40);m_ypout.precision(40);
        m_yout<< "y[" << x.second<<"]=" <<'\t' <<"("<< x.first[0] << ", " << x.first[1]<<"I)"<< "\n\n";
        m_ypout<<"y'[" << x.second<<"]=" <<'\t' <<"("<< x.first[2] << ", " << x.first[3]<<"I)"<< "\n\n";

        // }
      
  
   
    }
};


int main( int argc , char **argv )
{

   
    // Required precision of computations in ***decimal*** digits
    // Play with it to check different precisions
    const int digits = 50;

    // Setup default precision for all subsequent computations
    // MPFR accepts precision in bits - so we do the conversion
    mpreal::set_default_prec(mpfr::digits2bits(digits));
    
    
    //[initialisation and other params

    vector_type x(4);

    //use quotations to avoid loss of precison when initialising.
    
    //***************CHECK ICs MATCH INITIAL RADIUS*********************************
    x[0] = "-0.000089442307556012242046870383549907733735339434109588";
    x[1]= "0.000044722294418505882281368894833928366063067168385832" ;
    x[2] = "-4.4641753542932442508693361966916403862686093351869E-6";
    x[3] = "-8.9504832483903066707703454060478359939305628018336E-6";
    

    my_float dr = "-0.01"; //use neg time steps to go inwards
    my_float start = "10000.0";
    my_float end = "10";
    my_float omega = "0.1";
    my_float ell = "1.0";

    mpreal::set_default_prec(mpfr::digits2bits(40));  
    cout.precision(40); // Show only some of the digits

    my_float abs_err = "1.0E-40" , rel_err = "1.0E-25" , a_x = "1.0" , a_dxdt = "1.0";

    cout << abs_err << endl;
    cout << rel_err << endl;
    //]
     
 

    //[***********************rosenbrock (integrate working) ******************
     
    /*
      size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< my_float > >( abs_err , rel_err ) ,
      make_pair( radMod(omega , ell) , radMod_jacobi(omega , ell)) ,
      x , start , end , dr ,
      cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
    */
      
    //]

    //[***********************rosenbrock (integrate adaptive, working but not dense enough) ******************
     
    /*
      size_t num_of_steps = integrate_adaptive( make_dense_output< rosenbrock4< my_float > >( abs_err , rel_err ) ,
      make_pair( radMod(omega , ell) , radMod_jacobi(omega , ell)) ,
      x , start , end , dr ,
      cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
    */
      
    //]
      

    //[***********rb5 adaptive(works but won't accept 2 errs)*********
    
    
    typedef rosenbrock4< my_float > stepper_type;
    auto stepper = make_dense_output< stepper_type >( abs_err , rel_err );

    ofstream y,yp;
    y.open ("ydata.txt");
    yp.open ("ypdata.txt");
    y.precision(40);
    yp.precision(40);
    std::for_each(
        make_adaptive_time_iterator_begin( stepper ,
                                          make_pair( radMod(omega , ell) , radMod_jacobi(omega , ell)) ,
                                          x , start , end , dr) ,
        make_adaptive_time_iterator_end( stepper ,
                                        make_pair( radMod(omega , ell) , radMod_jacobi(omega , ell)) , x ) ,
        streaming_observer(y,yp)
        );

    


    return 0;
}