tode.cpp

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#include <iostream>
#include <iomanip>
#include <cmath>

#include "numerical/integschemesctes.h"
#include "numerical/autostepme.h"

class ODE
{
public:
    ODE()
    {
        x  = 0.0;
        y  = 0.5;
        Dx = 2.0;
        z  = 0;
    }
    int Solve()
    {
        IntegSchemesCtes isc; isc.ME_maxSS(4000).ME_STOL(1.0e-7);
        AutoStepME<double,double,double,ODE> TI(isc);
        return TI.Solve(this, &ODE::tangent_dy_dz, &ODE::calc_error,
                        x,y,z,Dx);
    }
    double x;
    double y;
    double Dx;
    double z;
private:
    int tangent_dy_dz(double const &  x, double const &  y, double const & z, 
                       double const & dx, double       & dy, double       & dz) const
    {
        // dy/dt=y-t^2+1 {0<=t<=2,y(0)=0.5} ==> x=t, z=0
        // solution: y(t)=(t+1)^2-0.5exp(t)
        dy = ( y-x*x+1 )*dx;
        dz = 0.0;
        return 0;
    }
    double calc_error(double const & Ey, double const & y_high,
                      double const & Ez, double const & z_high) const
    {
        return fabs(Ey)/fabs(y_high);
    }
}; // class ODE



using namespace std;

int main(int argc, char **argv)
{
    ODE ode;

    int k = ode.Solve(); // x will be x_1
    cout << "num of its (k) = " << k << endl;

    if (k==-1)
    {
        cout << "ERROR: TI did not converge!\n";
        return 1;
    }

    cout << "y_1 = " << setw(20) << setprecision(15) << ode.y << endl;
    cout << "solution:\n";
    cout << "y_1 = " << setw(20) << setprecision(15) << (ode.x+1.0)*(ode.x+1.0)-0.5*exp(ode.x) << endl;

    return 0;
}

---