Demo: Class ODESolver usage: A linear ODE given by numerical data

We present the same problem as the previous example but the differential equation is given as linear one by its coefficients. We consider the linear differential equation
 c1y'(t) + c0y(t) = f(t)
We test the code with the data:
 f(t)=2*exp(t)-1, c0=c1=1, y(0)=0.
With these data we obtain the solution y(t)=exp(t)-1.

  • The code stars in the same way as for the previous demo code.

    #include "OFELI.h"
    using namespace OFELI;
        
    int main(int argc, char *argv[])
    {
       theFinalTime = 1.;
       theTimeStep = atof(argv[1]);

  • The ODE is defined by declaring an instance of class ODESolver. The constructor defined the numerical scheme to solve the equation. Note that we use the enum variable TimeScheme that contains all implemented schemes. Other possiblities are FORWARD_EULER, FORWARD_EULER, HEUN and AB2. Note that only explicit methods are used. Implicit methods lead to a nonlinear equation at each time step, which we avoid here.

       ODESolver ode(RK4);

  • Once the instance created, we set the initial solution, and define then the ODE by the member function setF. We note that the used variables are t and y since we attempt to solve the equation
             y'(t) = F(t,y(t))
    It is important that the initial solution is to be given before defining the equation. We eventually run the time marching procedure.

       ODESolver ode(HEUN);
       ode.setInitial(0.);
       ode.setInitialRHS(1.);
                              
       TimeLoop {
          ode.setCoef(1,1,0,2*exp(theTime)-1);
          ode.runOneTimeStep();
       }
                              
       cout << ode << endl;
       cout << "Error: " << fabs(exp(theFinalTime)-1-ode.get()) << endl;
       return 0;
    }

  • We can now display the ode data and the solution. We also dispolay the error for this case where the exact solution is y(t) = exp(t)-1 .

    
       cout << ode << endl;
       cout << "Solution:  " << ode.get() << endl;
       cout << "Error:     " << fabs(exp(theFinalTime)-1-ode.get()) < endl;
       return 0;
    }