YANE-Framework Tutorial 1.1.0

examples/odesolve/src/solout_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 ofstream out;
00025 
00026 void dgl ( int * n, double * t, double * x, double * dx, double * rpar,
00027         int * ipar )
00028 {
00029         dx [ 0 ] = x [ 1 ];
00030         dx [ 1 ] = rpar [ 0 ] * ( 1 - x [ 0 ] * x [ 0 ] ) * x [ 1 ] - x [ 0 ];
00031 }
00032 
00033 void solout ( int * nr, double * t_old, double * t, double * x, int * n,
00034         double * cont, int * icont, int * nrd, double * rpar, int* ipar,
00035         int* ireturn )
00036 {
00037         out << * t << " " << x [ 1 ] << endl;
00038 }
00039 
00040 double rho ( int * n, double * t, double * y, double * rpar, int * ipar )
00041 {
00042         return ( rpar [ 0 ] + sqrt ( rpar [ 0 ] * rpar [ 0 ] - 4.0 ) ) / 2.0;
00043 }
00044 
00045 /***************************************************************************
00046  *                                                                         *
00047  * Main Routine                                                            *
00048  *                                                                         *
00049  ***************************************************************************/
00050 int main ( void )
00051 {
00052         double * rpar = new double [ 1 ];
00053         double * x0 = new double [ 2 ];
00054 
00055         //yane::OdeSolve::OdeConfig * odeconfig = new yane::OdeSolve::RodasConfig ( );
00056         //yane::OdeSolve::OdeSolveFirst * odesolver = new yane::OdeSolve::Rodas ( );
00057 
00058         //yane::OdeSolve::OdeConfig * odeconfig = new yane::OdeSolve::Radau5913Config ( );
00059         //yane::OdeSolve::OdeSolveFirst * odesolver = new yane::OdeSolve::Radau5913 ( );
00060 
00061         //yane::OdeSolve::OdeConfig * odeconfig = new yane::OdeSolve::DoPriConfig ( );
00062         //yane::OdeSolve::OdeSolveFirst * odesolver = new yane::OdeSolve::DoPri853 ( );
00063 
00064         yane::OdeSolve::OdeConfig * odeconfig = new yane::OdeSolve::OdeConfig ( );
00065         yane::OdeSolve::OdeSolveFirst * odesolver = new yane::OdeSolve::ROCK4 ( );
00066 
00067         yane::OdeSolve::OdeFunction * odefunc = new yane::OdeSolve::OdeFunction (
00068                 dgl, 2 );
00069 
00070         odefunc->setRho ( true, rho );
00071 
00072         out.open ( "data.txt" );
00073 
00074         cout << "Stiff example" << endl << endl;
00075         /*      rpar[0] = 5000.0;*/
00076         rpar [ 0 ] = 5.0;
00077         x0 [ 0 ] = 0.0;
00078         x0 [ 1 ] = 2.0;
00079 
00080         odefunc->setSolutionOutput ( solout );
00081 
00082         // Reset to current function and config
00083         odesolver->reset ( odefunc, odeconfig );
00084 
00085         // set t0, x0
00086         odesolver->init ( 0.0, x0 );
00087 
00088         // refresh config and solve the equation
00089         try
00090         {
00091                 odesolver->refreshConfig ( );
00092                 odesolver->calc ( 10.0, rpar );
00093                 cout << "Solution : x(" << 10 << ") = (" << x0 [ 0 ] << "," << x0 [ 1 ]
00094                         << ")" << endl;
00095         }
00096         catch ( yane::Utils::Exception e )
00097         {
00098                 cout << "Unable to solve equation." << endl;
00099         }
00100 
00101         cout << "Accepted steps: " << odesolver->acceptedSteps ( ) << endl;
00102 
00103         out.close ( );
00104 
00105         FILE * gnuplotpipe = popen ( "gnuplot", "w" );
00106 
00107         // Plot into file
00108         //fprintf(gnuplotpipe,"%s\n%s\n","set terminal png", "set output 'data.png'");
00109         fprintf ( gnuplotpipe, "%s\n", "plot  'data.txt'" );
00110         fflush ( gnuplotpipe );
00111 
00112         cin.get ( );
00113 
00114         fprintf ( gnuplotpipe, "exit\n" );
00115         pclose ( gnuplotpipe );
00116 
00117         delete odeconfig;
00118         delete odesolver;
00119         delete odefunc;
00120 
00121         delete [ ] rpar;
00122         delete [ ] x0;
00123 
00124         return 0;
00125 }