Ordinary differential equations (and initial value problems)

Ordinary differential equations (and initial value problems)#

In mathematics, an ordinary differential equation (ODE) is a differential equation (DE) dependent on only a single independent variable. As with other DE, its unknown(s) consists of one (or more) function(s) and involves the derivatives of those functions The term “ordinary” is used in contrast with partial differential equations which may be with respect to more than one independent variable.(see https://en.wikipedia.org/wiki/Ordinary_differential_equation?useskin=vector)

An initial value problem, IVP, is a a problem specified by an ordinary differential equation and its initial conditions.

You have been already solving a differential equation: The second Newton law

(44)#\[\begin{equation} ma = m\frac{dv}{dt} = m\frac{d^2x}{dt^2} = \sum F = f(t, x, v), \end{equation}\]

Can be seen as a first order differential equation (when talking about $v$) or as a second order differential equation (regarding $x$).

For instance, a projectile must follow a parabolic path, but that is because the parabolic path is the solution for the corresponding differential equation.

From: "https://upload.wikimedia.org/wikipedia/commons/thumb/a/a4/Parabolic_trajectory.svg/500px-Parabolic_trajectory.svg.png"

In this chapter we will focus on solving, numerically, Ordinary differential equations. Later we will learn some basic methods for Partial differential equations. For a nice introduction to both of them, please watch:

A numerical solution for and IVP corresponds to finding an algorithm that allows you to move from time $t$ to time $t+\delta t$, given a time discretization. It is assumed that the general problem is represented as

(45)#\[\begin{equation} \dot y = \frac{dy}{dt} = f(t, y), \end{equation}\]

where $y$ represents the state of the system. All linear differential equations of any order can be cast to this form by introducing auxiliary derivatives, so the problem becomes vectorial, but its form remains

(46)#\[\begin{equation} \frac{d\vec y}{dt} = \vec f(t, \vec y). \end{equation}\]

So the key question is how to go numerically from $\vec y(t)$ to $\vec y(t+\delta t)$. There are several algorithms, with both advantages and disadvantages (in sinplicity, precision, stability, etc), to do so. Here we will see Euler algorithm which is the most straight forward (but lacks precision and stability), and later one of the most common, the Runge-Kutta. After that, specific algorithms for second order equations will be discussed (like Leap-Frog and Verlet), and, finally some numerical libraries will be shown.

The Euler algorithm#

The Euler method is based on the forward derivative:

(47)#\[\begin{equation} \frac{d\vec y}{dt} = \frac{\vec y(t+h) - \vec y(t)}{h} + O(h). \end{equation}\]

Here, $h = \delta t$. Since we are looking for $f(t + \delta t)$, and we actually know the derivatives (they are given by the differential equation we want to solve), then

(48)#\[\begin{align} \vec y(t+\delta t) &= \vec y(t) + \delta t \frac{d \vec y}{dt} + O(\delta t^2), \\ &= \vec y(t) + \delta t \vec f(t, \vec y) + O(\delta t^2), \end{align}\]

Where $\vec f(t, \vec y) = d\vec y/dt$ is just the function returning the derivatives. As you can see , basically the Euler method is performing a linear extrapolation.

From: "https://upload.wikimedia.org/wikipedia/commons/thumb/1/10/Euler_method.svg/440px-Euler_method.svg.png"

Although easy to implement, it lacks stability, and its low order (locally $O(\delta t^2)$, globally $O(N\delta t^2) = O(\delta t)$) implies a very small $\delta t$, so many steps to solve the IVP on a given time interval. Check https://www.desmos.com/calculator/0i0oebdu6j

A simple, one dimensional implementation#

Here we will solve the following one dimensional differential equation using the Euler algorithm:

\[ \frac{dy}{dt} = f(y, \dot y, t) = r y (1-y), \]

where $r$ is a paremeter. This equation is called the “logistic equation” and used in the population dynamics field (https://en.wikipedia.org/wiki/Logistic_function?useskin=vector#Logistic_differential_equation). To solve this equation, we will need

a function that returns the derivatives: for now, it will just return a number. But how to handle the constant $r$? for now it will be just a global constant.
a function to perform the euler steps from and initial to a final time: It must also receive the initial state and then evolve it using the Euler algorithm.

A simple, inmediate implementation is s follows

#include <iostream>

const double R = 0.10;
double f(double y, double t);
void euler(double & y, double t0, double tf, double h);

int main(void)
{
  // condiciones iniciales
  double y = 0.57;

  // evolucion temporal
  euler(y, 0.0, 80.9876, 0.05);
  
  return 0;
}

double f(double y, double t)
{
  return R*y*(1-y);
}

void euler(double & y, double t0, double tf, double h)
{
  // avanzar desde y en t, hasta y en tf
  for(double t = t0; t <= tf; t = t + h) {
    std::cout << t << "\t" << y << "\t" << f(y, t) << "\n";
    y = y + h*f(y, t);
  }
}

Please compile it, run it and check the effect of $h$ and $r$.

A more general implementation: The Lorenz model and atractor#

The Lorenz model for climate represents a coupled system of one-dimensional differential equations that is used to study the climate and also gave the first evidence on sensitivy to initial conditions and the onset of chaos. Its solutions attractor is very famous:

From: "https://upload.wikimedia.org/wikipedia/commons/1/13/A_Trajectory_Through_Phase_Space_in_a_Lorenz_Attractor.gif"

It is represented as the following set of coupled equations:

(49)#\[\begin{align} \frac{dx}{dt} = \sigma(y-x),\\ \frac{dy}{dt} = x(\rho - z) - y, \\ \frac{dz}{dt} = xy - \beta z. \end{align}\]

As you can see, we have three independet variables, and three parameters.

Please about th possible modifications needed for the previous code to be able to solve the Lorens system. Some guiding questions:

How many dimensions are now in the state vector to be solved? how to represent that state?
how to return the derivative of a 3d vector?
How to pack the parameters to pass them to the derivatives and related functions?
How to print the system in a general way without mixing the integrator and the printing tools?

Exercise

Generalize the f function to be able to work on a N-dimensional state vector
Generalize the euler function to be able to work on a N-dimensional state vector
Pack the parameters in a given data structure.
Select some parameter values and plot two solutions that start very close.

General implementation of integration algorithms#

Before implementing the euler equation, and to easy future comparisons, let’s define the interface that we will use for all the implementations:

Main file

It must have the system constants, time interval and steps, etc.
System state type: state_t, will be represented by a data array, like a vector or valarray.
Derivatives function: It must return the derivatives, and must be defined in the main file since it corresponds to the actual differential equation you want to solve. But it could be vectorial, so we will include the ouput vector in the arguments. It models the function $\vec f(t, \vec y)$, so its definition will be void fderiv(const state_t & y, state_t & dydt, double t). Notice that the derivatives are a non-const reference, so we can store/return their values.

Header and implementation files: Functionality to solve an ODE The general prototype to solve the ODE will be a function with the following prototype

template <class deriv_t, class system_t, class printer_t>
    void integrate_euler(deriv_t deriv, system_t & s, double tinit, double tend, double dt, printer_t writer);

Why do we need to make it a template? because we want to allow the user to use any kind of type for the system class, and for printing the data. The first argument is the derivatives function, the second is the current state, then we have the nitial, final and time step, and finally something to write (wiether cout or a file stream).

This implementation might look complex but will help us to use other libs that use the same interface.

Please remember that the templates must go completely on the header file.

Euler algorithm implementation#

First of all, our Euler algorithm will be implemented as

template <class deriv_t, class system_t, class printer_t>
void integrate_euler(deriv_t deriv, system_t & y, double tinit, double tend, double dt, printer_t writer)
{
  // initial setup
  system_t dydt; // to store the derivatives
  const int N = y.size();
  dydt.resize(N);
  double time = tinit;
  int nsteps = (tend - tinit)/dt;
  writer(y, time); // write initial conditions
    
  // time evolution
  for(int ii = 1; ii < nsteps; ++ii) {
    time = tinit + ii*dt;
    deriv(y, dydt, time); // compute the new derivatives
    // y = y + dt*dydt; // coefficient wise (assumption)
    for (int ii = 0; ii < N; ++ii) { // for each component of the state vector
      y[ii] = y[ii] + dt*dydt[ii]; // euler algorithm
    }
    writer(y, time); // write the state for this new time
  }
}

Save this in the file integrator.h.

Solving a first order equation#

Please solve the first order differential equation

(50)#\[\begin{equation} \dot y = \frac{dy}{dt} = 5\sin(t) + t, \end{equation}\]

with initial condition $y(0) = 1.5$, for $t \in [0, 10.0]$. Plot the numerical and theoretical solution for $\delta t = 0.1, 0.2, 0.8, 1.5$. Our main file could be

#include <iostream>
#include <valarray>
#include <cmath>
#include <cstdlib>
#include "integrator.h"

typedef std::valarray<double> state_t; // alias for state type

void print(const state_t & y, double time);
void fderiv(const state_t & y, state_t & dydt, double t);

int main(int argc, char **argv)
{
  if (argc != 4) {
    std::cerr << "ERROR. Usage:\n" << argv[0] << " T0 TF DT" << std::endl;
    std::exit(1);
  }
  const double T0 = std::atof(argv[1]);
  const double TF = std::atof(argv[2]);
  const double DT = std::atof(argv[3]);

  const int N = 1;
  state_t y(N);

  // initial conditions
  y[0] = 1.5;

  // perform the actual integration
  // adapt to print to a file
  integrate_euler(fderiv, y, T0, TF, DT, print);

  return 0;
}

void fderiv(const state_t & y, state_t & dydt, double t)
{
    //
}

void print(const state_t & y, double time)
{
  std::cout << time << "\t" << y[0] << "\t" << y[1] << std::endl;
}

To compile just use the same command and everything will be fine

g++ -g -std=c++17 -fsanitize=address,undefined main_euler1.cpp -o main_euler1.x

Solving a second order differential equation#

Here we will simulate an harmonic oscillator, described by

(51)#\[\begin{equation} \frac{d^2x}{dt^2} + \omega^2 x = 0, \end{equation}\]

which can be transformed into a system of 2 differential equations using the velocity, $v = dx/dt$, as $\vec y = (y_0, y_1) = (x, dx/dt) = (x, v)$,

(52)#\[\begin{align} \frac{dy_0}{dt} &= \frac{dx}{dt} = v = y_1\\ \frac{dy_1}{dt} &= \frac{dv}{dt} = -\omega^2 x = -\omega^2 y_0. \end{align}\]

In this particular case we will need a parameter, $\omega$. In general, remember that if you need parameters (mass, frequency, spring constant) you either use lambda functions or functors to model the system derivatives and keep the derivative interface uniform. For instance, our derivative funtion could be a lambda like

double W = ...;
// actual derivatives: this is the system model
auto fderiv = [W](const state_t & y, state_t & dydt, double t){
    dydt[0] = y[1];
    dydt[1] = -W*W*y[0];
};

Solve the system, for $\omega = 1.234$ rad/s. Choose an appropriate time interval so you have at least for full oscillations. Again, plot the theoretical and numerical solutions for $\delta t = 0.1, 0.2, 0.5, 1.0$. Furthermore, plot the mechanical energy as function of time, and also plot the phase space $v$ versus $x$.

Exercises#

Use the previous code to check if the energy is conserved (plot the phase space, $v$ versus $x$). Check the effect of $\Delta t$.
Read, from a file, the parameters $\Delta t$, $t_0$, $t_f$, $\omega$.

Runge-Kutta algorithm#

Unfortunately, the Euler method is not stable, and its precision is far from ideal. We can improve its order by using methods that or order 3, 4 and so on. Here, we will study the typical Runge-Kutta method of order 4, for N coupled linear differential equations. This is much better than simple Euler. The price to pay is to evaluate the derivatives four times per time step.

From: "https://upload.wikimedia.org/wikipedia/commons/thumb/7/7e/Runge-Kutta_slopes.svg/600px-Runge-Kutta_slopes.svg.png"

For the RK4 algorithm , we can compute the evolution as

(53)#\[\begin{equation} \vec y(t+\delta t) = \vec y (t) + \frac{h}{6}(\vec k_1 + 2\vec k_2 + 2\vec k_3 + \vec k_4) + O(\delta t^5), \end{equation}\]

with the $\vec k_i$ functions evaluating the derivatives at different stages

(54)#\[\begin{align} \vec k_1 &= \vec f(t, \vec y(t)),\\ \vec k_2 &= \vec f(t + \frac{h}{2}, \vec y(t) + h\frac{\vec k_1}{2}),\\ \vec k_3 &= \vec f(t + \frac{h}{2}, \vec y(t) + h\frac{\vec k_2}{2}),\\ \vec k_4 &= \vec f(t + h, \vec y(t) + h\vec k_3),\\ \end{align}\]

The total global order is $O(\delta t^4)$.

According to the Runge-Kutta algorithm, we will need to go step by step computing k1, k2, ... and then update the state. We will compute the derivatives four times per time step, but the gain in precision will be huge, compared to the Euler.

To implement this, we just add a new function to our integrator.cpp file, adding the following implementation

template <class deriv_t, class system_t, class printer_t>
void integrate_rk4(deriv_t deriv, system_t & y, double tinit, double tend, double dt, printer_t writer)
{
    // init method
    system_t k1, k2, k3, k4, aux;
    const int N = y.size();
    k1.resize(N);
    k2.resize(N);
    k3.resize(N);
    k4.resize(N);
    aux.resize(N);
    // time evolution
    double time = tinit;
    int nsteps = (tend - tinit)/dt;
    writer(y, time); // write initial conditions
    for(int ii = 1; ii < nsteps; ++ii) {
        time = tinit + ii*dt;
        // k1 -> (t, y)
        deriv(y, k1, time);
        // k2 -> t+h/2, y+h*k1/2
        for (int ii = 0; ii < N; ++ii)
            aux[ii] = y[ii] + k1[ii]*dt/2;
        deriv(aux, k2, time + dt/2);
        // k3 -> t+h/2, y+h*k2/2
        for (int ii = 0; ii < N; ++ii)
            aux[ii] = y[ii] + k2[ii]*dt/2;
        deriv(aux, k3, time + dt/2);
        // k4 -> t+h, y+h*k3
        for (int ii = 0; ii < N; ++ii)
            aux[ii] = y[ii] + k3[ii]*dt;
        deriv(aux, k4, time + dt);
        // use all k
        for (int ii = 0; ii < N; ++ii)
            y[ii] = y[ii] + dt*(k1[ii] + 2*k2[ii] + 2*k3[ii] + k4[ii])/6.0;
        // call writer
        writer(y, time);
    }
}

Exercises#

Compare the results for the same system and same $\Delta t$ using both Euler and rk4. Does rk4 conserves the mechanical energy better?
Make the oscillator non-linear, with a force of the form $F(x) = -kx^\lambda$, with $\lambda$ even. How is the system behaviour as you change $\lambda$?
Add some damping to the linear oscillator, with the form $f_{damp} = -mbv$, where $m$ is the mass and $b$ the damping constant. Plot the phase space. Comment about the solution.
Solve many dynamical system that you can find in any numerical programming book.
Implement the leap-frog method: https://en.wikipedia.org/wiki/Leapfrog_integration?useskin=vector. Notice that you need a spetial treatment to start the algorithm.
Implement the Verlet method: https://en.wikipedia.org/wiki/Verlet_integration?useskin=vector . Notice that you need a spetial treatment to start the algorithm.

Using a library : odeint#

You could also use more elaborate libraries tailored to this kind of problems. In this section we will use odeint, https://headmyshoulder.github.io/odeint-v2/examples.html . Install it if you don’t have it. In the computer room it is already installed. The interface is pretty similar to what we have been doing, but it includes much more advanced algorithms (this is analogous to using scipy). Use the following example as a base code

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdlib>
#include <boost/numeric/odeint.hpp>

typedef std::vector<double> state_t; // alias for state type

void initial_conditions(state_t & y);
void print(const state_t & y, double time);
void fderiv(const state_t & y, state_t & dydt, double t);

int main(int argc, char **argv)
{
  if (argc != 5) {
    std::cerr << "ERROR. Usage:\n" << argv[0] << " T0 TF DT W" << std::endl;
    std::exit(1);
  }
  const double T0 = std::atof(argv[1]);
  const double TF = std::atof(argv[2]);
  const double DT = std::atof(argv[3]);
  const double W  = std::atof(argv[4]);
  const int N = 2;

  auto fderiv = [W](const state_t & y, state_t & dydt, double t){
    dydt[0] = y[1];
    dydt[1] = -W*W*y[0];
  };

  state_t y(N);
  initial_conditions(y);
  boost::numeric::odeint::integrate(fderiv, y, T0, TF, DT, print);

  return 0;
}

void initial_conditions(state_t & y)
{
  y[0] = 0.9876;
  y[1] = 0.0;
}

void print(const state_t & y, double time)
{
  std::cout << time << "\t" << y[0] << "\t" << y[1] << std::endl;
}

Exercises#

Again, compare the errors among euler, rk4 and odeint.
Use a functor to encapsulate the global parameter.
Solve the free fall of a body but with damping proportional to the squared speed.

More numerical libs#

Odeint: Odeint is a modern and fast library that provides a variety of ODE solvers, including Runge-Kutta methods, Adams-Bashforth methods, and backward difference formulas. It is also very flexible and can be used with a variety of different data structures and numerical types.
- http://scholarpedia.org/article/Odeint_library
- https://headmyshoulder.github.io/odeint-v2/
- https://github.com/headmyshoulder/odeint-v2
- https://www.boost.org/doc/libs/1_82_0/libs/numeric/odeint/doc/html/index.html
gsl, https://www.gnu.org/software/gsl/doc/html/ode-initval.html: The GNU Scientific Library (GSL) is a general-purpose library that contains a variety of mathematical and scientific functions. It also includes a number of ODE solvers, including Runge-Kutta methods, Adams-Bashforth methods, and backward difference formulas. Here is an example:

include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

// Harmonic oscillator function
int func (double t, const double y[], double f[], void *params) {
    double k = *(double *)params;
    f[0] = y[1];
    f[1] = -k * y[0];
    return GSL_SUCCESS;
}

int main() {
    double k = 1.0; // Spring constant
    gsl_odeiv2_system sys = {func, NULL, 2, &k};
    
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);
    int i;
    double t = 0.0, t1 = 10.0;
    double y[2] = { 1.0, 0.0 }; // Initial conditions

    for (i = 1; i <= 100; i++) {
        double ti = i * t1 / 100.0;
        int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

        if (status != GSL_SUCCESS) {
            std::cout << "Error: status = " << status << std::endl;
            break;
        }

        std::cout << t << " " << y[0] << " " << y[1] << std::endl;
    }

    gsl_odeiv2_driver_free (d);
    return 0;
}

sundials, https://computing.llnl.gov/projects/sundials: Sundials is another library for solving initial value problems for ordinary differential equations (ODEs). It is written in C and provides solvers for both non-stiff and stiff systems. Sundials includes CVODE, which is a solver for non-stiff problems, and IDA, which is a solver for stiff problems. Sundials also offers sensitivity analysis capabilities for systems with parametric dependencies.
- https://sundials.readthedocs.io/en/latest/
- https://computing.llnl.gov/projects/sundials/sundials-software Here is an examples for using this library (see also https://github.com/basilwong/simple-sundials/):

// 
    #include <iostream>
#include <sundials/sundials_types.h>
#include <sundials/sundials_math.h>
#include <sundials/sundials_nvector.h>
#include <sunmatrix/sunmatrix_dense.h> // access to dense SUNMatrix
#include <sunlinsol/sunlinsol_dense.h> // access to dense SUNLinearSolver
#include <cvode/cvode.h>
#include <nvector/nvector_serial.h>
#include <cvode/cvode_direct.h>

// Function to compute the right-hand side of the harmonic oscillator equation
int harmonicOscillator(realtype t, N_Vector y, N_Vector ydot, void *user_data) {
    realtype k = 1.0; // Spring constant
    NV_Ith_S(ydot, 0) = NV_Ith_S(y, 1); // y' = ydot
    NV_Ith_S(ydot, 1) = -k * NV_Ith_S(y, 0); // y'' = -k * y
    return 0;
}

int main() {
  int flag;
  const int N = 2;
  // Create the NVECTOR object
  N_Vector y = N_VNew_Serial(N); // 2 state variables

  // Set the initial conditions
  NV_Ith_S(y, 0) = 1.0; // Initial displacement
  NV_Ith_S(y, 1) = 0.0; // Initial velocity

  // Create the CVODE solver
  void *cvode_mem = CVodeCreate(CV_BDF);

  // Initialize the CVODE solver
  flag = CVodeInit(cvode_mem, harmonicOscillator, 0.0, y);
  std::clog << flag << "\n";
  
  // Set the integration parameters
  CVodeSStolerances(cvode_mem, 1e-6, 1e-8);
  CVodeSetMaxNumSteps(cvode_mem, 10000);

  // ---------------------------------------------------------------------------
  // Need to create a dense matrix for the dense solver.
  SUNMatrix A = SUNDenseMatrix(N, N);
  // ---------------------------------------------------------------------------

  // 9. Create Linear Solver Object.
  // ---------------------------------------------------------------------------
  // Dense linear solver object instead of the iterative one in the original
  // simple example.
  SUNLinearSolver LS = SUNDenseLinearSolver(y, A) ;
  // ---------------------------------------------------------------------------
  // ---------------------------------------------------------------------------
  // Call CVDlsSetLinearSolver to attach the matrix and linear solver this
  // function is different for direct solvers.
  flag = CVDlsSetLinearSolver(cvode_mem, LS, A);

  
  // Perform the integration
  realtype t_end = 10.0; // End time
  realtype t_out = 0.0;
  realtype t = 0.0;
  realtype dt = 0.1;
  for (t_out = dt; t_out <= t_end; t_out += dt) {
    CVode(cvode_mem, t_out, y, &t, CV_NORMAL);
    // Print the solution
    std::cout << t << "\t" << NV_Ith_S(y, 0) << "\t" << NV_Ith_S(y, 1) << std::endl;
  }
  // Free memory
  N_VDestroy(y);
  CVodeFree(&cvode_mem);

  return 0;
}

Eigen: Eigen is a C++ template library for linear algebra that can be used to solve initial value problems. It provides a variety of linear algebra operations and algorithms, including matrix operations, vector operations, and numerical solvers. Eigen’s numerical solvers can be used to solve systems of ordinary differential equations by converting them into matrix form and solving the resulting linear system[1].
NR, https://numerical.recipes/: Numerical Recipes is a popular book series that provides recipes for a variety of numerical algorithms. The book series also includes a number of ODE solvers, which are implemented in C++.
ODEPACK, https://computing.llnl.gov/projects/odepack: ODEPACK is a collection of Fortran solvers for ordinary differential equation (ODE) initial value problems. It includes solvers for both stiff and non-stiff systems of ODEs. LSODE is one of the solvers in ODEPACK that can handle both stiff and non-stiff systems. It uses Adams methods for non-stiff problems and Backward Differentiation Formula (BDF) methods for stiff problems. LSODE is available in separate double and single precision versions[3].
- https://computing.llnl.gov/projects/odepack

Modelling a bouncing ball and an introduction to OOP#

Now, we would like to model a ball bouncing on the floor. In this case, the derivatives function will get more complicated since the force with the floor is not present always, only when there is a collision with it. To model this interaction, we could follow the traditional way to do it by using the overlapping lenght $\delta$, which for a ball of radius $r$ and height $R_z$ can be computed as $\delta = r - R_z$. A collision with the wall will happen only when $\delta > 0$. What could be now the derivative function?

      // TODO: derivatives lambda

Now imagine that we want to simulate the same particle but in three dimensions. Then our state vector will need 6 components (three for the position, three for the velocity). What if we need 2, or 3, or $N$ particles? then we will have some disadvantages by using our current approach, like:

The state vector will need to have $6N$ components.
The actual indexing will be cumbersome.
The current force computation is done several times per time step (like in RK4, where we need 4 force evaluations). This will be slow for systems where computing the force is expensive.

To solve these problems in modelling, we will change our approach to more convenient Object Oriented Approach for the particles, and we will also use some integration algorithms better suited for this kind of problemas, like the leap-frog or verlet integration algorithms. This will be done in the next chapter

Ordinary differential equations (and initial value problems)

Contents

Ordinary differential equations (and initial value problems)#

The Euler algorithm#

A simple, one dimensional implementation#

A more general implementation: The Lorenz model and atractor#

General implementation of integration algorithms#

Euler algorithm implementation#

Solving a first order equation#

Solving a second order differential equation#

Exercises#

Runge-Kutta algorithm#

Exercises#

Using a library : odeint#

Exercises#

More numerical libs#

Modelling a bouncing ball and an introduction to OOP#