[tRock]

Hi,

I am trying to code the Lotka-Volterra equations ODE in OF and compare it with python https://docs.scipy.org/doc/scipy/ref...rate.solve_ivp

Based on the example in applications/test/ODE

I created:

#include "argList.H"
#include "IOmanip.H"
#include "ODESystem.H"
#include "ODESolver.H"
 
using namespace Foam;
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
class predatorPray
:
    public ODESystem
{
public:
    predatorPray(){} 
 
    label nEqns() const  {return 2;} 
 
    void derivatives
    (
        const scalar t,
        const scalarField& X,
        scalarField& dXdt
    ) const
    {
        const scalar a = 1.5;
        const scalar b = 1;
        const scalar c = 3;
        const scalar d = 1;
 
        const scalar x = X[0];
        const scalar y = X[1];
 
        dXdt[0] = a*x - b*x*y;
        dXdt[1] = -c*y  + d*x*y;
    }
 
    void jacobian
    (
        const scalar x,
        const scalarField& y,
        scalarField& dfdx,
        scalarSquareMatrix& dfdy
    ) const
    {
        FatalErrorInFunction << "Jacobian not defined " << abort(FatalError);
    }
};
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
 
int main(int argc, char *argv[])
{
 
    argList::addArgument("ODESolver");
    argList args(argc, argv);
 
    // Create the ODE system
    predatorPray ode;
 
    dictionary dict;
    dict.add("solver", args[1]);
 
    // Create the selected ODE system solver
    autoPtr<ODESolver> odeSolver = ODESolver::New(ode, dict);
 
    // Initial conditions
    scalarField dXdt(ode.nEqns());
    scalarField X(ode.nEqns());
    X[0] = 10;
    X[1] = 5;
 
    // Integration properties
    scalar tStart = 0;
    scalar tEnd = 15;
    scalar dt = 0.01;
 
    odeSolver->solve(tStart, tEnd, X, dt); // Default is relTol 1e-4 and aTol SMALL
 
    Info<< Foam::setprecision(12);
    Info << X << endl;
 
    Info<< "\nEnd\n" << endl;
 
    return 0;
}

I am using RKDP45 solver.


I would like to know if it is possible to check the integration value and current time value at each step of integration (So that I can also plot the data).


In addition, is it possible to evaluate the data at discrete time intervals as with python's solve_ivp?



Kind regards

[tRock]

Can anyone give me a hand with this? Would greatly appreciate it.