Examples ======== The following table gives an overview over all examples. ====================================================================================================== ================================================================== ================================================================================================================================= Section Brief Description File ====================================================================================================== ================================================================== ================================================================================================================================= `Arenstorf orbit <#arenstorf-orbit>`__ This example shows how to use multiple adaptive time step methods `arenstorf.cpp `__ `Brownian movement <#brownian-movement>`__ This example shows how to use random in an ODE `brownian.cpp `__ `Brusselator equations <#brusselator-equations>`__ The chemistry example of Brusselator (2 equations model) `brusselator.cpp `__ `Brusselator equations with DIRK method <#brusselator-equations-with-dirk-method>`__ This example shows hos to use DIRK methods `brusselator_dirk.cpp `__ `Curtiss-Hirschfelder equation <#curtiss-hirschfelder-equation>`__ This example shows how to use range and iterators on solution `curtiss_hirschfelder.cpp `__ `Curtiss-Hirschfelder equation with expRK method <#curtiss-hirschfelder-equation-with-exprk-method>`__ This example shows how to use exponential Runge-Kutta methods `curtiss_hirschfelder_exprk.cpp `__ `Exponential function <#exponential-function>`__ This example is the simplest example `exp.cpp `__ `Exponential function with exact solver <#exponential-function-with-exact-solver>`__ An example of splitting and exact solver `exp_splitting.cpp `__ `Heat model <#heat-model>`__ The classical heat equation solving with RKC2 method `heat.cpp `__ `ROCK method <#rock-method>`__ This example shows how to use ROCK2 and ROCK4 methods `heat_rock.cpp `__ `Samurai is hot <#samurai-is-hot>`__ This example shows how to coupling ponio and samurai `heat_samurai.cpp `__ `Lorenz equations <#lorenz-equations>`__ The chaotic system example of Lorenz equations `lorenz.cpp `__ `Lorenz equations with multiple methods <#lorenz-equations-with-multiple-methods>`__ This example shows how to use splitting methods and Lawson methods `lorenz_tuto.cpp `__ `Lorenz equations with all methods <#lorenz-equations-with-all-methods>`__ This example shows how to use all methods (except expRK) `lorenz_all.cpp `__ `Lotka-Volterra model <#lotka-volterra-model>`__ The classical predator–prey model of Lotka-Volterra `lotka_volterra.cpp `__ `Nagumo equation <#nagumo-equation>`__ Example of a traveling wave `nagumo.cpp `__ `Pendulum equation <#pendulum-equation>`__ The classical pendulum equation `pendulum.cpp `__ `Belousov Zhabotinsky <#belousov-zhabotinsky>`__ Solves Belousov-Zhabotinsky equations with PIROCK method `belousov_zhabotinsky_pirock.cpp `__ `Belousov Zhabotinsky in 2D <#belousov-zhabotinsky-in-2D>`__ Solves Belousov-Zhabotinsky equations with PIROCK method in 2D `bz_2d_pirock.cpp `__ ====================================================================================================== ================================================================== ================================================================================================================================= To lunch examples, in the main directory of ponio run: :: cmake . -B build -DBUILD_DEMOS=ON Eventually to get examples with samurai run: :: cmake . -B build -DBUILD_DEMOS=ON -DBUILD_SAMURAI_DEMOS=ON and :: cd build next you can compile an example with :: make AN_EXAMPLE you could also launch the Python script which launch the example and display results :: make AN_EXAMPLE_visu or :: make visu_AN_EXAMPLE Arenstorf orbit --------------- The system of differential equations for the Arenstorf orbit are: .. math:: \begin{cases} \ddot{x} &= x + 2\dot{y} - \frac{1-\mu}{r_1^3}(x+\mu) - \frac{\mu}{r_2^3}(x-1+\mu) \\ \ddot{y} &= y - 2\dot{x} - \frac{1-\mu}{r_1^3}y - \frac{\mu}{r_2^3}y \end{cases} where .. math:: r_1 = \sqrt{(x+\mu)^2 + y^2},\quad r_2 = \sqrt{(x-1+\mu)^2 + y^2} parameter :math:`\mu=0.012277471` and the initial condition gives by: .. math:: x(0) = 0.994,\quad \dot{x}(0) = 0,\quad y(0) = 0,\quad \dot{y}(0) = -2.00158510637908252240537862224 To solve this kind of problem with ponio, first of all you should rewrite it as the form: :math:`\dot{u} = f(t, u)`, here we classically take .. math:: u = \begin{pmatrix} x \\ y \\ \dot{x} \\ \dot{y} \end{pmatrix} So we have: .. math:: \dot{u} = \begin{pmatrix} \dot{x} \\ \dot{y} \\ \ddot{x} \\ \ddot{y} \end{pmatrix} = \begin{pmatrix} \dot{x} \\ \dot{y} \\ x + 2\dot{y} - \frac{1-\mu}{r_1^3}(x+\mu) - \frac{\mu}{r_2^3}(x-1+\mu) \\ y - 2\dot{x} - \frac{1-\mu}{r_1^3}y - \frac{\mu}{r_2^3}y \end{pmatrix} = f(t, u) In this example we solve this system with some explicit adaptive time step methods from `Dormand, J.R., Prince, P.J., A family of embedded Runge-Kutta formulae (1980) Journal of Computational and Applied Mathematics `__ ================= ==================== Arenstorf orbit Arenstorf velocity ================= ==================== |Arenstorf orbit| |Arenstorf velocity| ================= ==================== =================== ================= Time step history Error over time =================== ================= |Time step history| |Error over time| =================== ================= All example in `arenstorf.cpp `__, and run :: make arenstorf_visu Brownian movement ----------------- We write a simple Brownian movement .. math:: \begin{cases} \dot{x} = X(t) \\ \dot{y} = Y(t) \end{cases} where :math:`X(t)` and :math:`Y(t)` are random variable (juste a ``std::rand`` at each iteration). +------------------------------+ | Some Brownian movement in 2D | +==============================+ | |brownian movement| | +------------------------------+ All example in `brownian.cpp `__, and run :: make brownian_visu Brusselator equations --------------------- The Brusselator is a model of periodic chemical reaction. We present the version ODE with two species .. math:: \begin{cases} \dot{x} &= m_a - (m_b + 1)x + x^2y \\ \dot{y} &= m_b x - x^2y \end{cases} We solve this model with a hight order explicit Runge-Kutta method: RK(8, 6). =========================== ========================================== Brusselator concentration Brusselator concentration in phase space =========================== ========================================== |Brusselator concentration| |Brusselator concentration in phase space| =========================== ========================================== All example in `brusselator.cpp `__, and run :: make brusselator_visu Brusselator equations with DIRK method -------------------------------------- The Brusselator is a model of periodic chemical reaction. We present the version ODE with two species .. math:: \begin{cases} \dot{x} &= m_a - (m_b + 1)x + x^2y \\ \dot{y} &= m_b x - x^2y \end{cases} In this example we choose to solve the model with a diagonal-implicit Runge-Kutta method. The problem object has to be an ``ponio::implicit_problem`` and we need to compute the Jacobian matrix and proposes some linear algebra routines. For that we use `Eigen library `__. If ``state_t`` is floating point, a Eigen vector or a samurai field, ponio provides functions to solve implicit problems. In all cases, you can specify your own linear algebra object that contains a ``solver`` method (see ``lin_alg_2_2`` structure), that takes a matrix :math:`A` (same type as the returns type of jacobian gives to ``ponio::implicit_problem``) and a vector :math:`b` (same type as ``state_t``) and return the solution of the linear problem .. math:: Ax = b. ========================= ======================================== Brusselator concentration Brusselator concentration in phase space ========================= ======================================== |image1| |image2| ========================= ======================================== All example in `brusselator_dirk.cpp `__, and run :: make brusselator_dirk_visu Curtiss-Hirschfelder equation ----------------------------- A classical stiff problem is the Curtiss-Hirschfelder equation .. math:: \dot{y} = k(\cos(t) - y) with :math:`k>1` and :math:`y(0) = y_0`. We choose :math:`k = 50` and :math:`y_0 = 2`. In this example we present how to control time loop with a ``ponio::solver_range``. You can do it with an iterator on this range with : .. code:: cpp auto sol_range = ponio::make_solver_range( ... ); and iterate over this range with a classical iterator with: .. code:: cpp for ( auto it = sol_range.begin(); it < sol_range.end(); ++it ) { // ... // current time : it->time // current state : it->state // current time step : it->time_step } or with a range-based for loop: .. code:: cpp for ( auto ui : sol_range ) { // ... } Only in the first case you can control time step before increment (with your adaptive time step heuristic) with modification of ``it->time_step`` data member. +---------------------------------+ | Curtiss-Hirschfelder solution | +=================================+ | |Curtiss-Hirschfelder solution| | +---------------------------------+ All example in `curtiss_hirschfelder.cpp `__, and run :: make curtiss_hirschfelder_visu Curtiss-Hirschfelder equation with expRK method ----------------------------------------------- A classical stiff problem is the Curtiss-Hirschfelder equation .. math:: \dot{y} = k(\cos(t) - y) with :math:`k>1` and :math:`y(0) = y_0`. We choose :math:`k = 50` and :math:`y_0 = 2`. In this example we solve the equation with Krogstad method (an exponential Runge-Kutta method), and LRK(4, 4) method (a Lawson method). In both methods, you need to define a ``ponio::lawson_problem`` with a linear and non-linear part. We choose the linear part :math:`-k` and the non-linear part as :math:`N:t, y\mapsto k\cos(t)`. Exponential Runge-Kutta methods and Lawson methods are build to solve exactly the linear part when the non-linear part goes to 0. +-------------------------------+ | Curtiss-Hirschfelder solution | +===============================+ | |image3| | +-------------------------------+ All example in `curtiss_hirschfelder_exprk.cpp `__, and run :: make curtiss_hirschfelder_exprk_visu Exponential function -------------------- In this example we solve the simplest differential equation: .. math:: \dot{y} = y with :math:`y(0) = 1`. After a lot of calculus we can find the exact solution :math:`y(t) = e^t`, or we can approximate it with RK NSSP (2, 1) method. +------------------------+ | Exponential function | +========================+ | |Exponential function| | +------------------------+ All example in `exp.cpp `__, and run :: make exp_visu Exponential function with exact solver -------------------------------------- In this example we solve the simplest differential equation: .. math:: \dot{y} = \lambda y + (1-\lambda )y with :math:`y(0) = 1` and :math:`\lambda = 0.3`. In this example we solve this equation with: - a Strang splitting method with :math:`f_1:(t,y)\mapsto \lambda y` (which is solved exactly), and :math:`f_2:(t,y)\mapsto (1-\lambda)y` (which is solved with a RK(2,2) Raslton). - an exact solver :math:`\phi:(f,t,u,\Delta t)\mapsto e^{\Delta t}u`. - a RK(2,2) Ralston method. +----------------------+ | Exponential function | +======================+ | |image4| | +----------------------+ All example in `exp_splitting.cpp `__, and run :: make exp_splitting_visu Heat model ---------- In this example we propose to solve a PDE, the heat equation in 1D .. math:: \partial_t u = -\partial_{xx} u with the initial condition gives by the fundamental solution of head equation at time :math:`t=0.001` given by: .. math:: u(t, x) = \frac{1}{2\sqrt{\pi t}} e^{-\frac{x^2}{4t}} In ponio, the ``state_t`` should propose arithmetic operations as addition and multiplication by a scalar (of type ``value_t``). For the sake of simplicity, we use in the example a ``std::valarray``. The heat equation is quite complicated to solve with an explicit Runge-Kutta method but we do it with a extended stability method with the Runge-Kutta Chebyshev of order 2. In ponio you could choose the number of stages of this method : ``ponio::runge_kutta::explicit_rkc2<15>()`` (for 15 stages). +-----------------------------+ | Solution of heat equation | +=============================+ | |Solution of heat equation| | +-----------------------------+ All example in `heat.cpp `__, and run :: make heat_visu ROCK method ----------- In this example we propose to solve a PDE, the heat equation in 1D .. math:: \partial_t u = -\partial_{xx} u with the initial condition gives by the fundamental solution of head equation at time :math:`t=0.001` given by: .. math:: u(t, x) = \frac{1}{2\sqrt{\pi t}} e^{-\frac{x^2}{4t}} In ponio, the ``state_t`` should propose arithmetic operations as addition and multiplication by a scalar (of type ``value_t``). For the sake of simplicity, we use in the example a ``std::valarray``. An optimization of RKC2 is the ROCK2 method from `Abdulle, A., Medovikov, A. Second order Chebyshev methods based on orthogonal polynomials. Numer. Math (2001) `__, and its extension to order 4, ROCK4 method presented in `Abdulle, A. Fourth Order Chebyshev Methods with Recurrence Relation. SIAM Journal on Scientific Computing (2002) `__. ========================= ================================== Solution of heat equation Mesure of order of ROCK2 and ROCK4 ========================= ================================== |image5| |Mesure of order| ========================= ================================== All example in `heat_rock.cpp `__, and run :: make heat_rock_visu Samurai is hot -------------- This example needs to activate ``-DBUILD_SAMURAI_DEMOS=ON`` In this example we propose to solve a PDE, the heat equation in 1D .. math:: \partial_t u = -\partial_{xx} u with the initial condition gives by the fundamental solution of head equation at time :math:`t=0.001` given by: .. math:: u(t, x) = \frac{1}{2\sqrt{\pi t}} e^{-\frac{x^2}{4t}} In this example we coupling the mesh refinement library `samurai `__ with ponio. +--------------------------------------------------------+ | Solution of heat equation with levels of adaptive mesh | +========================================================+ | |image6| | +--------------------------------------------------------+ All example in `heat_samurai.cpp `__, and run :: make heat_samurai_visu Lorenz equations ---------------- The classical chaotic system example of Lorenz equations .. math:: \begin{cases} \dot{x} &= \sigma (y - x) \\ \dot{y} &= \rho x - y - xz \\ \dot{z} &= xy - \beta z \end{cases} In this example we use a ``std::vector`` to store the current state :math:`(x, y, z)` vector. ================ ======================= Solution in 3D Solution by composant ================ ======================= |Solution in 3D| |Solution by composant| ================ ======================= All example in `lorenz.cpp `__, and run :: make lorenz_visu Lorenz equations with multiple methods -------------------------------------- The classical chaotic system example of Lorenz equations .. math:: \begin{cases} \dot{x} &= \sigma (y - x) \\ \dot{y} &= \rho x - y - xz \\ \dot{z} &= xy - \beta z \end{cases} ============== ===================== Solution in 3D Solution by composant ============== ===================== |image7| |image8| ============== ===================== All example in `lorenz_tuto.cpp `__, and run :: make lorenz_tuto_visu Lorenz equations with all methods --------------------------------- The classical chaotic system example of Lorenz equations .. math:: \begin{cases} \dot{x} &= \sigma (y - x) \\ \dot{y} &= \rho x - y - xz \\ \dot{z} &= xy - \beta z \end{cases} +----------------+ | Solution in 3D | +================+ | |image9| | +----------------+ All example in `lorenz_all.cpp `__, and run :: make lorenz_all_visu Lotka-Volterra model -------------------- The classical predator–prey model of Lotka-Volterra: .. math:: \begin{cases} \dot{x} = \alpha x - \beta xy \\ \dot{y} = \delta xy - \gamma y \end{cases} with parameters :math:`\alpha=\frac{2}{3}`, :math:`\beta=\frac{4}{3}`, :math:`\gamma = \delta = 1`, and with the initial condition :math:`(x, y) = (x_0, x_0)`, with different values of :math:`x_0`. ======================= ========================= Prey predator history Solution in phase space ======================= ========================= |Prey predator history| |Solution in phase space| ======================= ========================= All example in `lotka_volterra.cpp `__, and run :: make lotka_volterra_visu Nagumo equation --------------- The Nagumo equation is a propagation of a traveling wave .. math:: \partial_t u = d \partial_{xx}u + ku^2(1-u) with parameter :math:`k=1`, :math:`d=1`, with the initial solution given by exact solution to time :math:`t=0`: .. math:: u(t, x) = \frac{\exp(c(x - x_0 - vt))}{1 + \exp(c(x-x_0 - vt))} where .. math:: v = \frac{1}{\sqrt{2}}\sqrt{kd},\quad c = -\frac{1}{\sqrt{2}}\sqrt{\frac{k}{d}},\quad x_0 = -10. ======================== ================ Solution with RKC(20, 2) Absolute error ======================== ================ |Nagumo solution| |Absolute error| ======================== ================ All example in `nagumo.cpp `__, and run :: make nagumo_visu Pendulum equation ----------------- The second order differential equation for the angle :math:`\theta` of a pendulum acted on by gravity with friction can be written: .. math:: \ddot{\theta} + b\dot{\theta} + c\sin(\theta) = 0 where :math:`b` and :math:`c` are positive constants, we take :math:`b=0.25`, :math:`c=5.0`. As `Arenstorf orbit <#arenstorf-orbit>`__, we have to rewrite problem as :math:`\dot{u} = f(t, u)` : .. math:: \partial_t \begin{pmatrix} \theta \\ \omega \end{pmatrix} = \begin{pmatrix} \omega \\ -b\omega - c\sin(\theta) \end{pmatrix} +--------------------------------------------+ | Pendulum equation (solved with RK (4,4)) | +============================================+ | |Pendulum equation (solved with RK (4,4))| | +--------------------------------------------+ All example in `pendulum.cpp `__, and run :: make pendulum_visu Belousov Zhabotinsky -------------------- This example needs to activate ``-DBUILD_SAMURAI_DEMOS=ON`` The Belousov Zhabotinsky reaction is a chemical periodic reaction, in this example we look at a 1D reduction model given by .. math:: \begin{cases} \partial_t a &= D_a \partial_{xx}a + \frac{1}{\mu}(-qa -ab + fc) \\ \partial_t b &= D_b \partial_{xx}b + \frac{1}{\epsilon}(qa - ab + b(1-b)) \\ \partial_t c &= D_c \partial_{xx}c + b -c \\ \end{cases} with parameter .. math:: \epsilon = 10^{-3},\quad \mu = 10^{-5},\quad f=3,\quad q = 2\cdot 10^{-4} and diffusion coefficients .. math:: D_a = \frac{1}{400},\quad D_b = \frac{1}{400},\quad D_c = \frac{0.6}{400} In this example we coupling the mesh refinement library `samurai `__ with ponio, and the system is solved by PIROCK method. +-------------------------------------------------------+ | Belousov Zhabotinsky system solved by PIROCK method | +=======================================================+ | |Belousov Zhabotinsky system solved by PIROCK method| | +-------------------------------------------------------+ All example in `belousov_zhabotinsky_pirock.cpp `__, and run :: make belousov_zhabotinsky_pirock_visu Belousov Zhabotinsky in 2D -------------------------- This example needs to activate ``-DBUILD_SAMURAI_DEMOS=ON`` The Belousov Zhabotinsky reaction is a chemical periodic reaction, in this example we look at a 1D reduction model given by .. math:: \begin{cases} \partial_t a &= D_a \Delta a + \frac{1}{\mu}(-qa -ab + fc) \\ \partial_t b &= D_b \Delta b + \frac{1}{\epsilon}(qa - ab + b(1-b)) \\ \partial_t c &= D_c \Delta c + b -c \\ \end{cases} with parameter .. math:: \epsilon = 10^{-2},\quad \mu = 10^{-5},\quad f=1.6,\quad q = 2\cdot 10^{-3} and diffusion coefficients .. math:: D_a = 2.5\cdot 10^{-3},\quad D_b = 2.5\cdot 10^{-3},\quad D_c = 1.5\cdot 10^{-3} In this example we coupling the mesh refinement library `samurai `__ with ponio, and the system is solved by PIROCK method. ================================ ==================================== Solution :math:`b` and :math:`c` Levels of adapted mesh on solution ================================ ==================================== |Solution b and c| |Levels of adapted mesh on solution| ================================ ==================================== All example in `bz_2d_pirock.cpp `__, and run :: make bz_2d_pirock_visu .. |Arenstorf orbit| image:: img/arenstorf/01.png .. |Arenstorf velocity| image:: img/arenstorf/02.png .. |Time step history| image:: img/arenstorf/03.png .. |Error over time| image:: img/arenstorf/04.png .. |brownian movement| image:: img/brownian/01.png .. |Brusselator concentration| image:: img/brusselator/01.png .. |Brusselator concentration in phase space| image:: img/brusselator/02.png .. |image1| image:: img/brusselator_dirk/01.png .. |image2| image:: img/brusselator_dirk/02.png .. |Curtiss-Hirschfelder solution| image:: img/curtiss_hirschfelder/01.png .. |image3| image:: img/curtiss_hirschfelder_exprk/01.png .. |Exponential function| image:: img/exp/01.png .. |image4| image:: img/exp_splitting/01.png .. |Solution of heat equation| image:: img/heat/01.png .. |image5| image:: img/heat_rock/01.png .. |Mesure of order| image:: img/heat_rock/02.png .. |image6| image:: img/heat_samurai/01.png .. |Solution in 3D| image:: img/lorenz/01.png .. |Solution by composant| image:: img/lorenz/02.png .. |image7| image:: img/lorenz_tuto/01.png .. |image8| image:: img/lorenz_tuto/02.png .. |image9| image:: img/lorenz_all/01.gif .. |Prey predator history| image:: img/lotka_volterra/01.png .. |Solution in phase space| image:: img/lotka_volterra/02.png .. |Nagumo solution| image:: img/nagumo/01.png .. |Absolute error| image:: img/nagumo/02.png .. |Pendulum equation (solved with RK (4,4))| image:: img/pendulum/01.png .. |Belousov Zhabotinsky system solved by PIROCK method| image:: img/belousov_zhabotinsky_pirock/01.png .. |Solution b and c| image:: img/bz_2d_pirock/01.png .. |Levels of adapted mesh on solution| image:: img/bz_2d_pirock/02.png