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

Computational Science Asked by kostas1335 on January 5, 2021

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_nright) + 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-hfrac{ 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-hfrac{ 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-hfrac{ 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 gright)(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 $$qleft( (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 = fleft( 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)$$;

Correct answer by Laurent90 on January 5, 2021