YANE-Framework Tutorial 1.1.0

examples/odesolve/src/benchmark_demo.cpp

00001 /***************************************************************************
00002  *                                                                         *
00003  * Copyright (C) 2011 by www.nonlinearmpc.com                              *
00004  *                                                                         *
00005  * Authors:                                                                *
00006  *  Juergen Pannek <juergen.pannek@googlemail.com>                         *
00007  *  Michael Schamel <michael.schamel@uni-bayreuth.de>                      *
00008  *                                                                         *
00009  ***************************************************************************/
00010 
00011 #include <yane/odesolve.h>
00012 #include <yane/utils.h>
00013 
00014 #include <iostream>
00015 #include <cmath>
00016 #include <iomanip>
00017 #include <string>
00018 #include <cstdlib>
00019 #include <cstring>
00020 #include <fstream>
00021 
00022 using namespace std;
00023 
00024 void dgl ( int * n, double * t, double * x, double * dx, double * rpar,
00025         int * ipar )
00026 {
00027         double hsqr = rpar [ 0 ] * rpar [ 0 ];
00028         int N = ipar [ 0 ];
00029 
00030         dx [ 0 ] = ( x [ 1 ] - 2.0 * x [ 0 ] ) / hsqr;
00031         for ( int i = 1; i < N - 1; i++ )
00032         {
00033                 dx [ i ] = ( x [ i - 1 ] - 2.0 * x [ i ] + x [ i + 1 ] ) / hsqr;
00034         }
00035         dx [ N - 1 ] = ( - 2.0 * x [ N - 1 ] + x [ N - 2 ] ) / hsqr;
00036 }
00037 
00038 double rho ( int * n, double * t, double * y, double * rpar, int * ipar )
00039 {
00040         return 4.0 / ( rpar [ 0 ] * rpar [ 0 ] );
00041 }
00042 
00043 void solve ( yane::OdeSolve::OdeConfig * odeconfig,
00044         yane::OdeSolve::OdeSolveFirst * odesolver,
00045         yane::OdeSolve::OdeFunction * odefunc, double * x, double * rpar,
00046         int * ipar )
00047 {
00048         yane::Utils::RTClock * clock = new yane::Utils::RTClock ( );
00049 
00050         // Reset to current function and config
00051         odesolver->reset ( odefunc, odeconfig );
00052 
00053         // set t0, x0
00054         odesolver->init ( 0.0, x );
00055 
00056         clock->reset ( );
00057 
00058         for ( int i = 1; i < 20; i++ )
00059         {
00060                 // refresh config and solve the equation
00061                 try
00062                 {
00063                         odesolver->refreshConfig ( );
00064                         odesolver->calc ( 0.01 * ( double ) i, rpar, ipar );
00065                 }
00066                 catch ( yane::Utils::Exception e )
00067                 {
00068                         cout << setw ( 9 ) << "exc" << " | ";
00069                         cout.flush ( );
00070                         return;
00071                 }
00072         }
00073         cout << setw ( 9 ) << clock->elapsedSeconds ( ) << " | ";
00074         cout.flush ( );
00075 }
00076 
00077 /***************************************************************************
00078  *                                                                         *
00079  * Main Routine                                                            *
00080  *                                                                         *
00081  ***************************************************************************/
00082 int main ( void )
00083 {
00084         yane::Utils::Exception::enableDebugMessage ( true );
00085 
00086         int dim_array [ ] =
00087         { 10, 25, 50, 75, 100, 250, 500, 750, 1000 };
00088         double tol_array [ ] =
00089         { 1E-3, 1E-6, 1E-9 };
00090 
00091         int dimension;
00092 
00093         int * ipar = new int [ 1 ];
00094         double * rpar = new double [ 1 ];
00095 
00096         yane::OdeSolve::OdeFunction * odefunc;
00097         double * x0;
00098         double * x;
00099 
00100         yane::OdeSolve::OdeConfig * dopriconfig =
00101                 new yane::OdeSolve::DoPriConfig ( );
00102         yane::OdeSolve::OdeSolveFirst * doprisolver =
00103                 new yane::OdeSolve::DoPri853 ( );
00104 
00105         yane::OdeSolve::OdeConfig * radauconfig =
00106                 new yane::OdeSolve::Radau5913Config ( );
00107         yane::OdeSolve::OdeSolveFirst * radausolver =
00108                 new yane::OdeSolve::Radau5913 ( );
00109 
00110         yane::OdeSolve::OdeConfig * odeconfig = new yane::OdeSolve::OdeConfig ( );
00111         yane::OdeSolve::OdeSolveFirst * rock2solver = new yane::OdeSolve::ROCK2 ( );
00112         yane::OdeSolve::OdeSolveFirst * rock4solver = new yane::OdeSolve::ROCK4 ( );
00113 
00114         yane::OdeSolve::OdeConfig * rodasconfig =
00115                 new yane::OdeSolve::RodasConfig ( );
00116         yane::OdeSolve::OdeSolveFirst * rodassolver = new yane::OdeSolve::Rodas ( );
00117 
00118         cout << fixed;
00119 
00120         for ( int tol = 0; tol < 3; tol++ )
00121         {
00122                 dopriconfig->setTolerance ( tol_array [ tol ], tol_array [ tol ] );
00123                 radauconfig->setTolerance ( tol_array [ tol ], tol_array [ tol ] );
00124                 rodasconfig->setTolerance ( tol_array [ tol ], tol_array [ tol ] );
00125                 odeconfig->setTolerance ( tol_array [ tol ], tol_array [ tol ] );
00126 
00127                 cout << "~~~ tolerance: " << scientific << tol_array [ tol ] << fixed
00128                         << " ~~~" << endl;
00129                 cout << setw ( 5 ) << "dim" << " | " << setw ( 9 ) << "Dopri853"
00130                         << " | " << setw ( 9 ) << "Radau5913" << " | " << setw ( 9 )
00131                         << "ROCK2" << " | " << setw ( 9 ) << "ROCK4" << " | " << setw (
00132                         9 ) << "Rodas" << " | " << setw ( 9 ) << "ROCK2 rho" << " | "
00133                         << setw ( 9 ) << "ROCK4 rho" << " | " << endl;
00134                 cout
00135                         << "-------------------------------------------------------------------------------"
00136                         << endl;
00137 
00138                 for ( int dim = 0; dim < 9; dim++ )
00139                 {
00140                         // Get dimension
00141                         dimension = dim_array [ dim ];
00142 
00143                         // Create function
00144                         odefunc = new yane::OdeSolve::OdeFunction ( dgl, dimension );
00145 
00146                         x = ( double * ) malloc ( sizeof(double) * dimension );
00147 
00148                         // Create x0
00149                         x0 = new double [ dimension ];
00150                         for ( int i = 0; i < dimension; i++ )
00151                         {
00152                                 x0 [ i ] = - pow ( ( double ) i / ( ( double ) dimension - 1.0 )
00153                                         - 0.5, 2 ) + 0.25;
00154                         }
00155 
00156                         ipar [ 0 ] = dimension;
00157                         rpar [ 0 ] = 1.0 / ( double ) ( dimension - 1 );
00158 
00159                         cout << setw ( 5 ) << dimension << " | ";
00160                         cout.flush ( );
00161 
00162                         // DoPri853
00163                         memcpy ( x, x0, sizeof(double) * dimension );
00164                         solve ( dopriconfig, doprisolver, odefunc, x, rpar, ipar );
00165 
00166                         // Radau5913
00167                         memcpy ( x, x0, sizeof(double) * dimension );
00168                         solve ( radauconfig, radausolver, odefunc, x, rpar, ipar );
00169 
00170                         // ROCK2
00171                         memcpy ( x, x0, sizeof(double) * dimension );
00172                         solve ( odeconfig, rock2solver, odefunc, x, rpar, ipar );
00173 
00174                         // ROCK4
00175                         memcpy ( x, x0, sizeof(double) * dimension );
00176                         solve ( odeconfig, rock4solver, odefunc, x, rpar, ipar );
00177 
00178                         // Rodas
00179                         memcpy ( x, x0, sizeof(double) * dimension );
00180                         solve ( rodasconfig, rodassolver, odefunc, x, rpar, ipar );
00181 
00182                         odefunc->setRho ( true, rho );
00183 
00184                         // ROCK2 with rho
00185                         memcpy ( x, x0, sizeof(double) * dimension );
00186                         solve ( odeconfig, rock2solver, odefunc, x, rpar, ipar );
00187 
00188                         // ROCK4 with rho
00189                         memcpy ( x, x0, sizeof(double) * dimension );
00190                         solve ( odeconfig, rock4solver, odefunc, x, rpar, ipar );
00191 
00192                         cout << endl;
00193 
00194                         free ( x );
00195                         delete [ ] x0;
00196                         delete odefunc;
00197                 }
00198                 cout << endl;
00199         }
00200 
00201         delete dopriconfig;
00202         delete doprisolver;
00203         delete radauconfig;
00204         delete radausolver;
00205         delete odeconfig;
00206         delete rock2solver;
00207         delete rock4solver;
00208         delete rodasconfig;
00209         delete rodassolver;
00210 
00211         delete [ ] rpar;
00212         delete [ ] ipar;
00213 
00214         return 0;
00215 }