# 1.7 Stiffness and Implicit Methods

## 1.7.4 Newton-Raphson: Implement Implicit Methods on Nonlinear Problems

When the ODE's are nonlinear, implicit methods require the solution of a nonlinear system of algebraic equations at each iteration. To see this, consider the use of the trapezoidal method for a nonlinear problem,

 $v^{n+1} = v^ n + \frac{1}{2}{\Delta t}\left[f(v^{n+1},t^{n+1})+f(v^ n,t^ n)\right].$ (1.127)

We can define the following residual vector for the trapezoidal method,

 $R(w) \equiv w - v^ n -\frac{1}{2}{\Delta t}\left[f(w,t^{n+1})+f(v^ n,t^ n)\right].$ (1.128)

Thus, $$v^{n+1}$$ for the trapezoidal method is given by the solution of,

 $R(v^{n+1}) = 0,$ (1.129)

which is a nonlinear algebraic system of equations for $$v^{n+1}$$.

One of the standard methods for solving a nonlinear system of algebraic equations is the Newton-Raphson method. It begins with an initial guess for $$v^{n+1}$$ and solves a linearized version of $$R=0$$ to find a correction to the initial guess for $$v^{n+1}$$. So, define the current guess for $$v^{n+1}$$ as $$w^ m$$ where $$m$$ indicates the sub-iteration in the Newton-Raphson method. Note, common useage is to call the iterations in the Newton-Raphson solution for $$v^{n+1}$$ sub-iterations since these are iterations which occur within every time iteration from $$n$$ to $$n+1$$. To find the correction, $$\Delta w$$, where

 $w^{m+1} = w^ m + \Delta w,$ (1.130)

we linearize and solve the nonlinear residual equation,

 $$\displaystyle R(w^{m+1})$$ $$\displaystyle =$$ $$\displaystyle 0, \nonumber$$ (1.131) $$\displaystyle R(w^ m + \Delta w)$$ $$\displaystyle =$$ $$\displaystyle 0, \nonumber$$ (1.132) $$\displaystyle R(w^ m) + \left.\frac{\partial R}{\partial w}\right|_{w^ m} \Delta w$$ $$\displaystyle =$$ $$\displaystyle 0, \nonumber$$ (1.133) $$\displaystyle \left.\frac{\partial R}{\partial w}\right|_{w^ m} \Delta w$$ $$\displaystyle =$$ $$\displaystyle -R(w^ m). \label{equ:NRsystem}$$ (1.134)

This last line is a linear system of equations for the correction since $${\partial R}/{\partial w}$$ is a $$d \times d$$ matrix when the original ODE's are a system of $$d$$ equations. For example, for the trapezoidal method,

 $\left.\frac{\partial R}{\partial w}\right|_{w^ m} = I - \frac{1}{2}{\Delta t}\left.\frac{\partial f}{\partial w}\right|_{w^ m}.$ (1.135)

Usually, the initial guess for $$v^{n+1}$$ is the previous iteration, i.e. $$w^0 = v^ n$$. So, the entire iteration from $$n$$ to $$n+1$$ has the following form,

1. Set initial guess: $$w^0 = v^ n$$ and $$m=0$$.

2. Calculate residual $$R(w^ m)$$ and linearization $${\partial R}/{\partial w}|_{w^ m}$$.

3. Solve Equation 1.134 for $$\Delta w$$.

4. Update $$w^{m+1} = w^ m + \Delta w$$.

5. Check if $$R(w^{m+1})$$ is small. If not, set $$m \leftarrow m+1$$ and repeat steps 1 through 4.