Squashed 'third_party/boostorg/odeint/' content from commit 6ff2719
Change-Id: If4892e29c1a5e6cf3a7aa51486a2725c251b0c7d
git-subtree-dir: third_party/boostorg/odeint
git-subtree-split: 6ff2719b6907b86596c3d43e88c1bcfdf29df560
diff --git a/examples/bulirsch_stoer.cpp b/examples/bulirsch_stoer.cpp
new file mode 100644
index 0000000..0b18f46
--- /dev/null
+++ b/examples/bulirsch_stoer.cpp
@@ -0,0 +1,100 @@
+/*
+ * bulirsch_stoer.cpp
+ *
+ * Copyright 2011-2013 Mario Mulansky
+ * Copyright 2011-2012 Karsten Ahnert
+ *
+ * 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 <iostream>
+#include <fstream>
+#define _USE_MATH_DEFINES
+#include <cmath>
+
+#include <boost/array.hpp>
+#include <boost/ref.hpp>
+
+#include <boost/numeric/odeint/config.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
+#include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::array< double , 1 > state_type;
+
+/*
+ * x' = ( - x*sin t + 2 tan x ) y
+ * with x( pi/6 ) = 2/sqrt(3) the analytic solution is 1/cos t
+ */
+
+void rhs( const state_type &x , state_type &dxdt , const double t )
+{
+ dxdt[0] = ( - x[0] * sin( t ) + 2.0 * tan( t ) ) * x[0];
+}
+
+void rhs2( const state_type &x , state_type &dxdt , const double t )
+{
+ dxdt[0] = sin(t);
+}
+
+
+ofstream out;
+
+void write_out( const state_type &x , const double t )
+{
+ out << t << '\t' << x[0] << endl;
+}
+
+int main()
+{
+ bulirsch_stoer_dense_out< state_type > stepper( 1E-8 , 0.0 , 0.0 , 0.0 );
+ bulirsch_stoer< state_type > stepper2( 1E-8 , 0.0 , 0.0 , 0.0 );
+
+ state_type x = {{ 2.0 / sqrt(3.0) }};
+
+ double t = M_PI/6.0;
+ //double t = 0.0;
+ double dt = 0.01;
+ double t_end = M_PI/2.0 - 0.1;
+ //double t_end = 100.0;
+
+ out.open( "bs.dat" );
+ out.precision(16);
+ integrate_const( stepper , rhs , x , t , t_end , dt , write_out );
+ out.close();
+
+ x[0] = 2.0 / sqrt(3.0);
+
+ out.open( "bs2.dat" );
+ out.precision(16);
+ integrate_adaptive( stepper , rhs , x , t , t_end , dt , write_out );
+ out.close();
+
+ x[0] = 2.0 / sqrt(3.0);
+
+ out.open( "bs3.dat" );
+ out.precision(16);
+ integrate_adaptive( stepper2 , rhs , x , t , t_end , dt , write_out );
+ out.close();
+
+
+ typedef runge_kutta_dopri5< state_type > 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 = make_dense_output( 1E-9 , 1E-9 , dopri5_type() );
+
+ x[0] = 2.0 / sqrt(3.0);
+
+ out.open( "bs4.dat" );
+ out.precision(16);
+ integrate_adaptive( dopri5 , rhs , x , t , t_end , dt , write_out );
+ out.close();
+
+}