Backward differentiation formula

The Backward differentiation formula (BDF) is a family of implicit methods for the numerical integration of ordinary differential equations. They are linear multistep methods that, for a given function and time, approximate the derivate of that function using informationfrom already computed time points, thereby increasing the accuracy of the approximation.

1. General formuma

We use BDF to solve the initial value problem

\[\begin{equation*} y'=f(t,y), y(t_0)=y_0 \end{equation*}\]

The general formula for a BDF can be written as

\[\begin{equation*} \sum_{k=0}^s a_k y_{n+k} = h\beta f(t_{n+s},y_{n+s}) \end{equation*}\]

where h denotes he step size and \(t_n=t_0+nh\). Since \(f\) is evaluated for the unknown \(y_{n+s}\), BDF methods are implicit and possibly require the solution of nonlinear equations at each step. The coefficients \(\alpha_k\) and \(\beta\) are chosen so that the method achieves order s, which is the maximum possible.

2. Derivation of the ccoefficients

Starting from the formula \(y'(t_{n+s})=f(t_{n+s},y(t_{n+s}))\) one approximates \(y(t_{n+s})=y_{n+s}\) and \(y'(t_{n+s}) \approx p'_{n,s}(t_{n+s})\) where \(p'_{n,s}(t)\) is the Lagrange interpolation polynomial for the points \((t_n,y_n),\ldots,(t_{n+s},y_{n+s})\). Using that \(t_n=t_0+nh\) and multiplying by h one arrives at the BDF method of order s.

For exemple the BDF with order 1 is :

\[\begin{equation*} y_{n+1} - y_n = hf(t_{n+1},y_{n+1}) \end{equation*}\]

and order 2 is :

\[\begin{equation*} y_{n+2} - \frac{4}{3}y_{n+1}+\frac{1}{3}y_n = \frac{2}{3}hf(t_{n+2},y_{n+2}) \end{equation*}\]

3. BDF in feelpp

To implement BDF on feelpp, we use the feel bdf library.

First, we define our bdf, one for each space :

auto mybdfA = bdf(_space = Ah, _name="mybdfA");
auto mybdfV = bdf(_space = Vh, _name="mybdfV");

We must then initialize each bdf with the initial conditions :

  for (auto time : mybdfA -> priorTimes() )
  {
    if (Environment::worldComm().isMasterRank())
    {
      std::cout << "Initialize prior times (from timeInitial()) : " << time.second << "s index: " << time.first << "\n";
    }
    A0.setParameterValues({{"t",time.second}});
    A0e = project(_space=Ah, _expr=A0);
    mybdfA->setUnknown(time.first,A0e);
  }
  for (auto time : mybdfV -> priorTimes() )
  {
    if (Environment::worldComm().isMasterRank())
    {
      std::cout << "Initialize prior times (from timeInitial()) : " << time.second << "s index: " << time.first << "\n";
    }
    V0.setParameterValues({{"t",time.second}});
    V0e = project(_space=Vh, _expr=V0);
    mybdfV->setUnknown(time.first,V0e);
  }

We finally start our bdf right before the time loop.

  mybdfA->start();
  mybdfV->start();
  for (double t = dt; mybdfA->isFinished() == false; )
  {
      ...

We define in the loop the polynomial derivate coefficients associeted to mybdfA, because only the derivate of A appear in both equations.

auto bdfA_poly = mybdfA->polyDeriv();

Thereafter, the only changes in the loop compared to an explicit Euler formulation are at the level of the equations, and at the level of the time variable. Here is the difference between the two versions for the implementation of the first equation.

For explicit Euler it was :

	  M00 += integrate( _range=markedelements(mesh, material.meshMarkers()),
			    _expr = dt * 1/mur * trace(trans(gradt(A))*grad(A)) );

	  M00 += integrate( _range=markedelements(mesh, material.meshMarkers()),
			    _expr = mu0 * sigma * inner(id(A) , idt(A) ));

	  M01  += integrate(_range=markedelements(mesh, material.meshMarkers()),
			    _expr = dt * mu0 * sigma * inner(id(A),trans(gradt(V))) );

	  F0 += integrate(_range=markedelements(mesh, material.meshMarkers()),
			  _expr = mu0 * sigma * inner(id(A) , idv(Aold)));

and for bdf it’s:

	    M00 += integrate( _range=markedelements(mesh, material.meshMarkers()),
			      _expr = 1/mur * trace(trans(gradt(A))*grad(A)) );

	    M00 += integrate( _range=markedelements(mesh, material.meshMarkers()),
			      _expr = mu0 * mybdfA->polyDerivCoefficient(0) * sigma * inner(id(A) , idt(A) ));

	    M01  += integrate(_range=markedelements(mesh, material.meshMarkers()),
			      _expr = mu0 * sigma * inner(id(A),trans(gradt(V))) );

	    F0 += integrate(_range=markedelements(mesh, material.meshMarkers()),
			    _expr = mu0 * sigma * inner(id(A) , idv(bdfA_poly)));

Finally, instead of using \(t\) we use

mybdfA->time()

or

mybdfB->time()

to get the time.

At the end of the loop, we initialized next step of both bdf with :

    mybdfA->next(*A);
    mybdfV->next(*V);

The default order is 1, if you want to do the simulation with a higher order, the option --bdf.order n, where n is the order, can be added before execution. You can also add the field "order = n", in the bdf part of the cfg file associated with the configuration.

4. Results

To run a simulation, you can compile the mqs-form.cpp file in the bdf branch and use for example the command line :

mpirun -np 4 feelpp_mqs_form --config-file cases/quart-turn/quart-turn.cfg --gmsh.hsize=0.05 --pc-type gasm --ksp-monitor=1 --bdf.order 2

If we run simulations on the two test cases we saw earlier, we get the same results. The implementation of the bdf schema is therefore correct.