How to derive the simplified Newton iteration in the TR-BDF2 ODE integration scheme

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.

Bibliography

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.

• 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? Nov 17 '20 at 18:58
• 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. Nov 17 '20 at 19:07
• Working on the increment also requires only solving the Newton equations to a lower level of accuracy. Nov 17 '20 at 21:21
• Yes, $f$ is assumed sufficiently smooth and $h$ sufficiently small. Nov 18 '20 at 9:22
• @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 ? Nov 22 '20 at 15:11