BDF2: How to Implement

Question

While I had a course a couple years ago concerning the implementation of linear time marching methods, I'm pretty rusty and am getting stumped here. I'm attempting to implement the BDF2 time marching method in order to convect Lagrangian points in space and time. What I have are the locations of these points and the velocities at which the points need to be convected (which is really just linear advection), but the underlying methodology that governs the velocity at each point is said to be unstable for explicit TMMs and BDF2 was specifically recommended as a higher order method to ensure stability after the first time step.

So, the equation that I'm trying to integrate (where r is a position vector, t is time, and V is velocity):

$$ \frac{\delta \vec{r}}{\delta t}=\vec{V} $$

The RHS of this expression is highly nonlinear, although it is analytical and I have discrete values of V readily available.

So, the query: how do I implement a time marching method (specifically BDF2) for this problem? The formulation that I've found is as follows (where n denotes the time step level): $$ \vec{r}_{n+1}=\frac{4}{3}\vec{r}_n-\frac{1}{3}\vec{r}_{n-1}+\frac{2}{3}(\Delta t)\frac{\delta \vec{r}_{n+1}}{\delta t} $$

The kicker for me is discretizing the last term (i.e., the derivative of r with respect to t). My thought is that I could just make a 2nd order Taylor approximation for the derivative using a predicted value the future time step (using the Euler explicit method) and then the actual values at present and prior time steps and use that. Would that be a correct methodology?

Answer 1

1. General framework

I assume your PDE is of the form $$ \frac{\partial r}{\partial t}(t,x)=f(t,x,\mathcal{D}r(t,x)) $$ where $$ \mathcal{D}r=(r,Dr,D^{2}r) $$ and $Dr$ is the gradient of $r$ and $D^{2}r$ is the Hessian of $r$. Then, BDF2 is just $$ r_{i}^{n+1}-\frac{4}{3}r_{i}^{n}+\frac{1}{3}r_{i}^{n}=\frac{2}{3}f(t^{n+1},x_{i},(\mathcal{D}_{h}r^{n+1})_{i})h.\tag{*} $$ In the above, $h=\Delta t$ is the timestep size and $$ (\mathcal{D}_{h}r^{n+1})_{i}=(r_{i}^{n+1},(D_{h}r^{n+1})_{i},(D_{h}^{2}r^{n+1})_{i}) $$ where $D_{h}^{p}$ is the matrix containing the stencil for the $p$-th derivative. (*) is a nonlinear equation in $$\vec{r}^{n+1}=(r_{1}^{n+1},\ldots,r_{M}^{n+1})$$ and hence must be solved iteratively (e.g., by Newton's method).

Remark: While I have only included up to second-order derivatives, you can add more.

2. Concrete Example

Consider the PDE $$ \frac{\partial r}{\partial t}(t,x)=\frac{\partial^2 r}{\partial x^2}(t,x) $$ This is just a special case of the PDE in Section 1 (i.e., take $f(t,x,\mathcal{D}r(t,x))=\frac{\partial^2 r}{\partial x^2}(t,x)$). In this case, BDF2 looks like $$ r_{i}^{n+1}-\frac{4}{3}r_{i}^{n}+\frac{1}{3}r_{i}^{n}=\frac{2}{3}\frac{r^{n+1}_{i-1}-2r^{n+1}_i+r^{n+1}_{i+1}}{(\Delta x)^2}h $$ where we have used a central difference to discretize the second derivative.