The Problem

The TR-BDF2 explained in this paper [1], is quite a popular numerical scheme used to integrate $\dot{y} = f(t,y)$, consistent of the following two stages:

\begin{align} y_{n+\gamma} & = y_n + \gamma \frac{h}{2}\left( f_n + f_{n+\gamma} \right) \tag{1} \\ % y_{n+1} & = \frac{1}{\gamma(2-\gamma)}y_{n+\gamma} - \frac{(1-\gamma)^2}{\gamma(2-\gamma)}y_n + \frac{1-\gamma}{2-\gamma}hf_{n+1} \tag{2} \end{align}

The above are two implicit algebraic equations that need to be solved in each step. This is accomplished by performing a simplified Newton (chord) iteration, where in accordance to Hosea & Shampine's paper [1] is performed by firstly letting $z=hf(t,y)$. As a result by substituting this into equations $(1)$ and $(2)$ above for the $k^{\text{th}}$ Newton iteration, we have:

\begin{align} y_{n+\gamma}^k & = \left(y_n + \frac{\gamma}{2}z_n\right) + \frac{ \gamma} {2}z_{n+\gamma}^k \tag{3}\\ % y_{n+1}^k & = \left(y_n +\frac{\sqrt{2}}{4}z_n + \frac{\sqrt{2}}{4}z_{n+\gamma} \right) + \frac{\gamma}{2}z_{n+1}^k \tag{4} \end{align}

From the above, the paper goes straight away to say that the Newton step and Newton iteration are as follows:

\begin{align} \left(I-h\frac{ \gamma}{2}\frac{\partial f}{\partial y}\right)\Delta^k = hf(t_{n+\gamma},y^k_{n+\gamma}) - z_{n+\gamma}^k, \ \text{and} \ z_{n+\gamma}^{k+1} = z_{n+\gamma}^k + \Delta^k \ \text{for equation (3)} \tag{5} \\ % \left(I-h\frac{ \gamma}{2}\frac{\partial f}{\partial y}\right)\Delta^k = hf(t_{n+1},y^k_{n+1}) - z_{n+1}^k, \ \text{and} \ z_{n+1}^{k+1} = z_{n+1}^k + \Delta^k \ \text{for equation (4)} \tag{6} \end{align}

The Question

Therefore I do not understand how the last two equations above were derived, and what the quantity $z$ represents.

To elaborate, the paper says that $z=hf(t,y)$, so $z_{n+\gamma}^k=hf(t_{n+\gamma},y_{n+\gamma}^k)$, but doesn't that make the RHS of equations $(5)$ and $(6)$ zero? Also, the typical Newton iteration is of the form $u^{k+1} = u^k - (g')^{-1}g(u^k)$. In this case then what is the function $g$ corresponding to for equations $(5)$ and $(6)$? It surely can't be $g = hf(t_{n+\gamma},y^k_{n+\gamma}) - z_{n+\gamma}^k$, since the Jacobian of this function is not $I-h\frac{ \gamma}{2}\frac{\partial f}{\partial y}$ as shown in equation $(5)$ for example.

Thanks for your help in advance.


  1. Hosea, M. E.; Shampine, L. F., Analysis and implementation of TR-BDF2, Appl. Numer. Math. 20, No. 1-2, 21-37 (1996). ZBL0859.65076.

Maybe it would be easier if you first consider the Newton method on the variable $y_{n+\gamma}$. I'll only treat the first stage of the method which defines $y_{n+\gamma}$, as the second one can be handled similarly. I'll also remove the dependence on time as it superfluous for the explanation. Let's recall this first stage:

$$y_{n+\gamma} = y_n + \gamma \frac{h}{2}\left( f_n + f_{n+\gamma} \right)$$

Let's introduce $y=y_{n+\gamma}$ and reformulate as: $$ 0 = g(y) = y - \frac{\gamma h}{2} f(y) - \underbrace{\left(y_n + \frac{\gamma h}{2} f(y_n) \right)}_{A}$$ where $A$ is constant.

To solve this, we apply the Newton method on $g$, starting from an initial guess $y^0$ for $y$, and the $k$-th Newton step reads: $$y^{k} = y^{k-1} - J(y_k)^{-1} g(y^k)$$

with the Jacobian $J(y^k)= \left(\partial_y g\right)(y^k)$.

Replacing $J(y^k)$ by $J(y^0)$ yields the corresponding simplified Newton method. Otherwise the full Newton method would compute the Jacobian of $g$ at each iterate anew.

Now the authors of the reference you gave have introduced $z=f(y)$. The first stage can then be reformulated as: $$y = y_n + \gamma \frac{h}{2}\left( z_n + z \right)$$

That is not what we will solve, as it would not make sense on its own. Indeed we have introduced a new variable, hence an additional equation is required to solve the system. That is why the system to be solved is now: $$y = y_n + \gamma \frac{h}{2}\left( z_n + z \right)$$ $$z = f(y)$$ which is equivalent to $q\left( (y,z)^t \right)=0$

You could write the Newton method to solve $q=0$. But here $y$ can be explicitly expressed from $z$ from the first equation. Hence we take the RHS of the first equation and use it in the second to obtain: $$z = f\left( y_n + \gamma \frac{h}{2}\left( z_n + z \right)\right)$$ $$\Leftrightarrow 0 = g_2(z)$$

You apply a Newton loop on $g_2$ as shown before for $g$, but this time the Jacobian is $\partial_z g_2 = \mathrm{I} - \frac{\gamma h}{2} \partial_y f$. Thus we obtain Equations $(5)$ and $(6)$ from your post.

So, regarding your last question: you actually precisely want the RHS of $(5)$ and $(6)$ to got to 0 as your Newton method converges. Then, when convergence is achieved, you have $z=f(y)$, from which you can compute $y$ based on Equation $(5)$.

EDIT: for future reference, the article by Hosea and Shampine explains one advantage of solving for $z_{n+1}=(d_t y)(t_{n+1})$ instead of $y_{n+1}$ itself. As the final stage of one step of the method is the same as the first of the next step (FSAL = First Same As Last), if we solved for $y_{n+1}$, then we would, by recomputing $z$ as $f(y)$ "excite" the fast modes of the solution and lose a bit of the stability of the method. Therefore it is best to compute $z_{n+1}$ and then recompute $y_{n+1}$ from the Runge-Kutta quadrature. Then at the next step, we directly set $z_n$ equal to the previous $z_{n+1}$. This is well explained in the paper.

  • $\begingroup$ Thanks for your answer. Just to clarify $$g_2(z) = z - f\left( y_n + \gamma \frac{h}{2}\left( z_n + z \right)\right)$$, is that correct? $\endgroup$ – kostas1335 Nov 17 '20 at 18:58
  • $\begingroup$ Yes it is ! In other papers, you might find that the Newton method is performed not on $y$ or $z=f(y)$, but rather on $\Delta y = y - y_n$. This may also be more convenient to implement. $\endgroup$ – Laurent90 Nov 17 '20 at 19:07
  • $\begingroup$ Working on the increment also requires only solving the Newton equations to a lower level of accuracy. $\endgroup$ – Wolfgang Bangerth Nov 17 '20 at 21:21
  • $\begingroup$ Yes, $f$ is assumed sufficiently smooth and $h$ sufficiently small. $\endgroup$ – kostas1335 Nov 18 '20 at 9:22
  • $\begingroup$ @WolfgangBangerth I don't understand why the increment allows for a looser tolerance. Let's take the first stage of this BDF method. We can write the residual on the stage value $y$ as $g_1(y)=y-y_n-ahf(y)-hQ$, with $Q$ the term involving the previous step. If we write a residual on the increment $\delta=y-y_n$ instead, we get $g_2(\delta)=\delta- ahf(y_n+\delta)-hQ$. The Newton Jacobians $d_y g_1$ and $d_\delta g_2$ are the same, and the Newton steps will also be the same. The only difference I see is if we scale their norms differently for the stopping criterion. Am I missing something ? $\endgroup$ – Laurent90 Nov 22 '20 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.