YANE-Framework Tutorial 1.1.0

examples/odesolve/src/radaudopri_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         dx [ 0 ] = x [ 1 ];
00028         dx [ 1 ] = rpar [ 0 ] * ( 1 - x [ 0 ] * x [ 0 ] ) * x [ 1 ] - x [ 0 ];
00029 }
00030 
00031 void jacobi ( int * n, double * t, double * x, double * dfx, int * ldfx,
00032         double * rpar, int* ipar )
00033 {
00034         dfx [ 0 ] = 0;
00035         dfx [ 1 ] = 1;
00036         dfx [ 2 ] = 2 * rpar [ 0 ] * x [ 0 ] - 1;
00037         dfx [ 3 ] = rpar [ 0 ] * ( 1 - x [ 0 ] * x [ 0 ] );
00038 }
00039 
00040 void solve ( yane::OdeSolve::OdeConfig * odeconfig,
00041         yane::OdeSolve::OdeSolveFirst * odesolver,
00042         yane::OdeSolve::OdeFunction * odefunc, double * x, double * rpar )
00043 {
00044         // Reset to current function and config
00045         odesolver->reset ( odefunc, odeconfig );
00046 
00047         // set t0, x0
00048         odesolver->init ( 0.0, x );
00049 
00050         cout << setw ( 7 ) << "t" << " | " << setw ( 20 ) << "x" << endl;
00051         cout << "--------------------------------" << fixed << setprecision ( 5 )
00052                 << endl;
00053 
00054         for ( int i = 1; i < 20; i++ )
00055         {
00056                 // refresh config and solve the equation
00057                 try
00058                 {
00059                         odesolver->refreshConfig ( );
00060                         odesolver->calc ( 0.2 * ( double ) i, rpar );
00061                         cout << 0.2 * ( double ) i << " | ( " << setw ( 8 ) << x [ 0 ]
00062                                 << ", " << setw ( 8 ) << x [ 1 ] << " )" << endl;
00063                 }
00064                 catch ( yane::Utils::Exception e )
00065                 {
00066                         cout << "Unable to solve equation." << endl;
00067                 }
00068         }
00069         cout << "--------------------------------" << endl;
00070 }
00071 
00072 void compare ( double * x0, double * rpar )
00073 {
00074         // DoPri
00075         yane::Utils::RTClock * clock = new yane::Utils::RTClock ( );
00076         clock->reset ( );
00077 
00078         yane::OdeSolve::OdeConfig * odeconfig = new yane::OdeSolve::DoPriConfig ( );
00079         yane::OdeSolve::OdeSolveFirst * odesolver = new yane::OdeSolve::DoPri5 ( );
00080         yane::OdeSolve::OdeFunction * odefunc = new yane::OdeSolve::OdeFunction (
00081                 dgl, 2 );
00082 
00083         double * tmp_x = new double [ 2 ];
00084 
00085         tmp_x [ 0 ] = x0 [ 0 ];
00086         tmp_x [ 1 ] = x0 [ 1 ];
00087         solve ( odeconfig, odesolver, odefunc, tmp_x, rpar );
00088 
00089         cout << "runtime using DoPri: " << clock->elapsedSeconds ( ) << endl
00090                 << endl;
00091 
00092         delete odeconfig;
00093         delete odesolver;
00094         delete odefunc;
00095 
00096         // Radau
00097         clock->reset ( );
00098 
00099         odeconfig = new yane::OdeSolve::RadauConfig ( );
00100         odesolver = new yane::OdeSolve::Radau5 ( );
00101         odefunc = new yane::OdeSolve::OdeFunction ( dgl, 2 );
00102         odefunc->setJacobiMatrix ( jacobi );
00103 
00104         tmp_x [ 0 ] = x0 [ 0 ];
00105         tmp_x [ 1 ] = x0 [ 1 ];
00106         solve ( odeconfig, odesolver, odefunc, tmp_x, rpar );
00107 
00108         cout << "runtime using Radau: " << clock->elapsedSeconds ( ) << endl
00109                 << endl;
00110 
00111         delete odeconfig;
00112         delete odesolver;
00113         delete odefunc;
00114 
00115         delete clock;
00116 }
00117 
00118 /***************************************************************************
00119  *                                                                         *
00120  * Main Routine                                                            *
00121  *                                                                         *
00122  ***************************************************************************/
00123 int main ( void )
00124 {
00125         double * rpar = new double [ 1 ];
00126         double * x0 = new double [ 2 ];
00127 
00128         rpar [ 0 ] = 1.0;
00129         x0 [ 0 ] = 2.0;
00130         x0 [ 1 ] = 0.0;
00131 
00132         cout << "Minimal non-stiff example" << endl << endl;
00133         compare ( x0, rpar );
00134 
00135         cout << "Stiff example" << endl << endl;
00136         rpar [ 0 ] = 5000.0;
00137         x0 [ 0 ] = 2.0;
00138         x0 [ 1 ] = 0.0;
00139 
00140         compare ( x0, rpar );
00141 
00142         cout << "Another Stiff example" << endl << endl;
00143         rpar [ 0 ] = 10000.0;
00144         x0 [ 0 ] = 2.0;
00145         x0 [ 1 ] = 0.0;
00146 
00147         compare ( x0, rpar );
00148 
00149         delete [ ] rpar;
00150         delete [ ] x0;
00151 
00152         return 0;
00153 }
00154