Lotka-Volterra system with some exponential Runge-Kutta method
==============================================================

In this example we present how to write a program to solve the
Lotka-Volterra equations with Ponio.

The system is definied as:

.. math::


     \begin{aligned}
       \frac{\mathrm{d}x}{\mathrm{d}t} &= \alpha x  - \beta xy \\
       \frac{\mathrm{d}y}{\mathrm{d}t} &= \delta xy - \gamma y \\
     \end{aligned}

where:

- :math:`x` is the number of prey
- :math:`y` is the number of predators
- :math:`t` represents time
- :math:`\alpha`, :math:`\beta`, :math:`\gamma` and :math:`\delta` are
  postive real parameters describing the interaction of the two species.

We would like to compute the invariant :math:`V` definied by:

.. math::


     V = \delta x - \gamma \ln(x) + \beta y - \alpha \ln(y)

We will use the same decomposition into linear and non-linear as the
Lawson presentation:

.. math::


       \dot{u} = Lu + N(t,u)

with :math:`u=(x,y)` and:

.. math::


       L = \begin{pmatrix}\alpha & 0 \\ 0 & -\gamma\end{pmatrix}
       \qquad
       N:(t,u)\mapsto \begin{pmatrix}-\beta u_1u_2 \\ \delta u_1u_2 \end{pmatrix}

Since the linear part is diagonal, we can use a ``std::valarray`` to
reprensent it.

.. code:: ipython3

    %system mkdir -p lotka_volterra_expRK_demo




.. parsed-literal::

    []



.. code:: ipython3

    %%writefile lotka_volterra_expRK_demo/main.cpp
    
    #include <iostream>
    #include <valarray>
    
    #include "ponio/solver.hpp"
    #include "ponio/observer.hpp"
    #include "ponio/problem.hpp"
    #include "ponio/runge_kutta.hpp"
    
    int main()
    {
        using state_t = std::valarray<double>;
        using namespace ponio::observer;
    
        double alpha = 2./3., beta=4./3., gamma=1., delta=1.;
        
        state_t L = { alpha, -gamma };
        auto N = [=]( double t, state_t const& u, state_t& du ) {
            double x=u[0], y=u[1];
            du[0] = -beta*x*y;
            du[1] = delta*x*y;
        };
        auto pb  = ponio::make_lawson_problem(L,N);
    
        double dt = 0.1;
        double tf = 100;
        state_t u_ini = {1.8,1.8};
        
        ponio::solve(pb, ponio::runge_kutta::exprk22_t<double, state_t>()            , u_ini, {0.,tf}, dt, "lotka_volterra_expRK_demo/exprk_exprk22.dat"_fobs);
        ponio::solve(pb, ponio::runge_kutta::krogstad_t<double, state_t>()           , u_ini, {0.,tf}, dt, "lotka_volterra_expRK_demo/exprk_k.dat"_fobs);
        ponio::solve(pb, ponio::runge_kutta::cox_matthews_t<double, state_t>()       , u_ini, {0.,tf}, dt, "lotka_volterra_expRK_demo/exprk_cm.dat"_fobs);
        ponio::solve(pb, ponio::runge_kutta::strehmel_weiner_t<double, state_t>()    , u_ini, {0.,tf}, dt, "lotka_volterra_expRK_demo/exprk_sw.dat"_fobs);
        ponio::solve(pb, ponio::runge_kutta::hochbruck_ostermann_t<double, state_t>(), u_ini, {0.,tf}, dt, "lotka_volterra_expRK_demo/exprk_ho.dat"_fobs);
        ponio::solve(pb, ponio::runge_kutta::etd2cf3_t<double, state_t>()            , u_ini, {0.,tf}, dt, "lotka_volterra_expRK_demo/exprk_etd2cf3.dat"_fobs);
        ponio::solve(pb, ponio::runge_kutta::etd3rk_t<double, state_t>()             , u_ini, {0.,tf}, dt, "lotka_volterra_expRK_demo/exprk_etd3rk.dat"_fobs);
    
        return 0;
    }


.. parsed-literal::

    Writing lotka_volterra_expRK_demo/main.cpp


.. code:: ipython3

    %system $CXX -std=c++20 -I ../include lotka_volterra_expRK_demo/main.cpp -o lotka_volterra_expRK_demo/main




.. parsed-literal::

    []



.. code:: ipython3

    %system ./lotka_volterra_expRK_demo/main




.. parsed-literal::

    []



.. code:: ipython3

    import numpy as np
    import matplotlib.pyplot as plt
    
    methods = [
        ("exprk22", "expRK(2,2)"),
        ("k", "Krogstad"),
        ("cm", "Cox-Matthews"),
        ("sw", "Strehmel Weiner"),
        ("ho", "Hochbruck-Ostermann"),
        ("etd2cf3", "ETD2CF3"),
        ("etd3rk", "ETD3RK")
    ]
    
    
    for tag, name in methods:
        data = np.loadtxt(f"lotka_volterra_expRK_demo/exprk_{tag}.dat")
        x_e = data[:,1]
        y_e = data[:,2]
        plt.plot(x_e[:int(len(x_e)/10)], y_e[:int(len(x_e)/10)], "-", label=name)
        del data
    
    plt.title("Solution in phase plane")
    plt.legend()
    plt.show()



.. image:: lotka_volterra_expRK_files/lotka_volterra_expRK_7_0.png


.. code:: ipython3

    def V(x,y, alpha=2./3., beta=4./3., gamma=1., delta=1.):
        return delta*x - np.log(x) + beta*y - alpha*np.log(y)

.. code:: ipython3

    for i,(tag, name) in enumerate(methods):
        data = np.loadtxt(f"lotka_volterra_expRK_demo/exprk_{tag}.dat")
        t_e = data[:,0]
        x_e = data[:,1]
        y_e = data[:,2]
        plt.plot(t_e, V(x_e, y_e), "-", label=name)
        del data
    plt.title("time evolution of $V$")
    plt.legend()
    plt.show()



.. image:: lotka_volterra_expRK_files/lotka_volterra_expRK_9_0.png


The same graph fro only high order methods

.. code:: ipython3

    methods = [
        ("k", "Krogstad"),
        ("cm", "Cox-Matthews"),
        ("sw", "Strehmel Weiner"),
        ("ho", "Hochbruck-Ostermann")
    ]
    for i,(tag, name) in enumerate(methods):
        data = np.loadtxt(f"lotka_volterra_expRK_demo/exprk_{tag}.dat")
        t_e = data[:,0]
        x_e = data[:,1]
        y_e = data[:,2]
        plt.plot(t_e, V(x_e, y_e), "-", label=name)
        del data
    plt.title("time evolution of $V$")
    plt.legend()
    plt.show()



.. image:: lotka_volterra_expRK_files/lotka_volterra_expRK_11_0.png