Subject: [Boost-commit] svn:boost r57616 - sandbox/odeint/libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-11-12 15:51:19
Author: mariomulansky
Date: 2009-11-12 15:51:18 EST (Thu, 12 Nov 2009)
New Revision: 57616
URL: http://svn.boost.org/trac/boost/changeset/57616
Log:
added harmonic_oscillator.cpp to examples
Added:
sandbox/odeint/libs/numeric/odeint/examples/harmonic_oscillator.cpp (contents, props changed)
Added: sandbox/odeint/libs/numeric/odeint/examples/harmonic_oscillator.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/examples/harmonic_oscillator.cpp 2009-11-12 15:51:18 EST (Thu, 12 Nov 2009)
@@ -0,0 +1,292 @@
+/* Boost numeric/odeint/examples/harm_osc.cpp
+
+ Copyright 2009 A. Constantine Ashford
+
+ Demonstrates the use of odeint to estimate the (known) solution to a simple,
+ undriven harmonic oscillator:
+
+ y'' = -y - gamma*y'
+
+ This example uses the Euler solver as well as he fourth-order Runge Kutta.
+ The output of each is printed in the output. The accuracy can be compared, or
+ the program can be profiled for a speed comparison.
+
+ The analytical solution to this ODE is (with initial displacement y and no
+ initial velocity y'):
+
+ Y = exp( -.5*gamma*t ) * cos( sqrt(1-gamma^2)*t )
+
+ This example calculates the analytical solution using the standard math
+ library and outputs the values at the same time points for comparison and
+ plotting.
+
+ On a unix system, a plot of the output could be created easily using gnuplot
+ and the following commands:
+
+ $ plot "harm_osc.out" using 1:2 title 'Euler' with dots,\
+ $ "harm_osc.out" using 1:5 title 'Runge Kutta 4' with dots,\
+ $ "harm_osc.out" using 1:8 title 'Exact' with dots
+
+ The step size, number of steps to integrate, output file name, and gamma can
+ be specified on the command line IFF you have the boost library files installed
+ *which are at the same version as this example is compiled against*. You can
+ enable the command line options by compiling with the preprocessor variable
+ BOOST_PROGRAM_OPTIONS_INSTALLED defined.
+
+
+ Compile with:
+ g++ -Wall -O3 -march=k8 -I$BOOST_DIR -I../../../../ harmonic_oscillator.cpp
+
+ 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>
+#include <iomanip>
+#include <vector>
+#include <cmath>
+
+#include <boost/numeric/odeint.hpp>
+
+#ifdef BOOST_PROGRAM_OPTIONS_INSTALLED
+#include <boost/program_options.hpp>
+namespace opts=boost::program_options;
+#endif
+
+namespace odeint=boost::numeric::odeint;
+
+// Namespace for variables specific to this example
+namespace harm_osc
+{
+ double gamma = .15;
+ double dt = 0.01;
+ size_t olen = 10000;
+ std::ofstream fout;
+
+ // The type of container used to hold the state vector
+ typedef std::vector<double> state_type;
+
+ /*
+ void harmonic_oscillator(state_type &x, state_type &dxdt, double t)
+
+ Calculates the derivatives of the state vector for the harmonic oscillator
+ system. The system equation is described in the opening comments.
+
+ The state vector is:
+ x[0] = Y (the output)
+ x[1] = Y' (dY/dt)
+
+ Arguments:
+ state_type &x: reference to the current state vector [ Y Y' ]
+ state_type &dxdt: reference to the vector of calculated derivitives
+ [ Y' Y'']
+ double t: the current integration time
+ */
+ void harmonic_oscillator(const state_type &x, state_type &dxdt, const double t)
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - gamma*x[1];
+ }
+
+ /*
+ Optionally included function to parse the command line. See detailed
+ comments and implementation at the bottom of this file
+ */
+ int parse_command_line(int ac, char ** av);
+}
+
+// Main
+int main(int argc, char **argv)
+{
+ using namespace std;
+
+ // Save the original destination of cout in case we redirect to a file
+ // as a result of --filename command line option
+ streambuf * cout_buf = cout.rdbuf();
+
+ // Parse the command line options iff boost::program_options is installed
+ #ifdef BOOST_PROGRAM_OPTIONS_INSTALLED
+ int parse_ret = harm_osc::parse_command_line(argc, argv);
+ if(parse_ret)
+ return parse_ret;
+ #endif
+
+ // Declare state vectors for the Euler system, and the Runge Kutta system
+ harm_osc::state_type x_euler(2);
+ harm_osc::state_type x_rk4(2);
+
+ // Set identical initial conditions for the Euler solver and the Runge Kutta
+ x_euler[0] = 1.0;
+ x_euler[1] = 0.0;
+ x_rk4[0] = 1.0;
+ x_rk4[1] = 0.0;
+
+ // Initialize the solvers
+ odeint::stepper_euler<harm_osc::state_type> euler;
+ odeint::stepper_rk4<harm_osc::state_type> rk4;
+
+
+ // Write the output header. The '#' symbol creates a comment in gnuplot
+ cout << "# Output from the simple harmonic oscillator example." << endl;
+ cout << "#" << endl;
+ cout << "# Columns:" << endl;
+ cout << "#" << endl;
+ cout << "# 1: time (s)" << endl;
+ cout << "# 2: Euler output" << endl;
+ cout << "# 3: Euler output dY/dt" << endl;
+ cout << "# 4: Euler output error" << endl;
+ cout << "# 5: 4th Order Runge Kutta output" << endl;
+ cout << "# 6: 4th Order Runge Kutta dY/dt" << endl;
+ cout << "# 7: 4th Order Runge Kutta output error" << endl;
+ cout << "# 8: Analytical (calculated) solution" <<endl;
+ cout << "#" << endl;
+ cout << "# A. Constantine Ashford" << endl;
+ cout << "#" << endl;
+
+// Initialize the time variable, and the analytical solution variable
+ double t = 0.0;
+ double reference;
+ double e_error = 0, r_error = 0;
+ double e_error_max = 0, r_error_max = 0;
+ double e_error_rms = 0, r_error_rms = 0;
+ int data_column_width = 15;
+
+ // Integrate and write the approximate and analytical results to a table
+ for(size_t oi=0; oi < harm_osc::olen; ++oi, t += harm_osc::dt)
+ {
+ // Calculate the analytical solution for this time
+ reference = exp(-.5*harm_osc::gamma*t) *
+ cos(sqrt(1-pow(harm_osc::gamma,2))*t);
+
+ // Calculate absolute errors for the data table
+ e_error = abs((reference - x_euler[0])/reference);
+ r_error = abs((reference - x_rk4[0])/reference);
+
+ // Write all the results to a row in the table
+ cout << setw(5) << t
+ << setw(data_column_width) << x_euler[0]
+ << setw(data_column_width) << x_euler[1]
+ << setw(data_column_width) << e_error
+ << setw(data_column_width) << x_rk4[0]
+ << setw(data_column_width) << x_rk4[1]
+ << setw(data_column_width) << r_error
+ << setw(data_column_width) << reference
+ << endl;
+
+ // Aggregate statistics
+ if( (e_error = abs(e_error)) > e_error_max)
+ e_error_max = e_error;
+
+ if((r_error = abs(r_error)) > r_error_max)
+ r_error_max = r_error;
+
+ e_error_rms += pow(e_error, 2);
+ r_error_rms += pow(r_error, 2);
+
+ // Advance the solvers
+ euler.next_step(harm_osc::harmonic_oscillator,
+ x_euler, t, harm_osc::dt);
+
+ rk4.next_step(harm_osc::harmonic_oscillator,
+ x_rk4, t, harm_osc::dt);
+
+ }
+
+ e_error_rms = sqrt(e_error_rms/harm_osc::olen);
+ r_error_rms = sqrt(r_error_rms/harm_osc::olen);
+
+ // Restore cout to direct to standard out in case we redirected to a file
+ cout.rdbuf(cout_buf);
+ cout << "#Integration complete." << endl;
+ cout << "#Errors for the Euler solver: MAX: " << e_error_max << " RMS: "
+ << e_error_rms << endl;
+ cout << "#Errors for the Runge Kutta solver: MAX: " << r_error_max << " RMS: "
+ << r_error_rms << endl;
+
+
+ // Success
+ return 0;
+}
+
+#ifdef BOOST_PROGRAM_OPTIONS_INSTALLED
+/*
+parse_command_line(int ac, char ** av)
+
+Uses boost::program_options to parse the command line.
+Allows the example to use runtime options such as:
+--step=.001 or
+--step .001
+
+Arguments:
+ int ac: the number of command line arguments. Usually just pass in argc
+ from main.
+
+ char ** av: an array of c-style strings holding the command line arguments.
+ Usually just pass in argv from main.
+
+Return value: 1 if the user chose the --help option, and does *not* wish to
+run an integration. Otherwise 0.
+
+*/
+int harm_osc::parse_command_line(int ac, char ** av)
+{
+ using namespace std;
+
+ // Declare supported options
+ opts::options_description opts_desc("Allowed options");
+ opts_desc.add_options()
+ ("help", "Display this message.")
+ ("file", opts::value<string>(), "Output file for data.(Default: stdout)")
+ ("step", opts::value<double>(), "Set step size (in seconds)")
+ ("duration", opts::value<size_t>(), "Set integration duration (in steps)")
+ ("gamma", opts::value<double>(), "Set gamma equation parameter");
+
+ // Let Boost Program Options parse the command line
+ opts::variables_map opts_map;
+ opts::store(opts::parse_command_line(ac, av, opts_desc), opts_map);
+ opts::notify(opts_map);
+
+ // Print help message then exit
+ if(opts_map.count("help"))
+ {
+ cout << opts_desc << endl;
+ return 1;
+ }
+
+ // Alter step size in seconds
+ if(opts_map.count("step"))
+ {
+ harm_osc::dt = opts_map["step"].as<double>();
+ cout << "Step size set to " << harm_osc::dt << " seconds" << endl;
+ }
+
+ // Change integration duration (in steps)
+ if(opts_map.count("duration"))
+ {
+ harm_osc::olen = opts_map["duration"].as<size_t>();
+ cout << "Integration duration set to " << harm_osc::olen << " steps."
+ << endl;
+ }
+
+ // Change damping constant
+ if(opts_map.count("gamma"))
+ {
+ harm_osc::gamma = opts_map["gamma"].as<double>();
+ cout << "Parameter GAMMA set to " << harm_osc::gamma << "." << endl;
+ }
+
+ // Redirect output from standard out to file
+ if(opts_map.count("file"))
+ {
+ cout << "Writing output to file " << opts_map["file"].as<string>()
+ << endl;
+ fout.open(opts_map["file"].as<string>().c_str());
+ cout.rdbuf(fout.rdbuf());
+ }
+
+ return 0;
+}
+#endif
+