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

Question

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.

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.

Answer 1

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.