Systems of nonlinear algebraic equations

Implicit time discretization methods for a system of ODEs, or a PDE, lead to systems of nonlinear algebraic equations, written compactly as

$$F(u) = 0,$$

where $u$ is a vector of unknowns $u=(u_0,\ldots,u_N)$, and $F$ is a vector function: $F=(F_0,\ldots,F_N)$. Sometimes the equation system has a special structure because of the underlying problem, e.g.,

$$A(u)u = b(u),$$

with $A(u)$ as an $(N+1)\times (N+1)$ matrix function of $u$ and $b$ as a vector function: $b=(b_0,\ldots,b_N)$.

We shall next explain how Picard iteration and Newton’s method can be applied to systems like $F(u)=0$ and $A(u)u=b(u)$. The exposition has a focus on ideas and practical computations. More theoretical considerations, including quite general results on convergence properties of these methods, can be found in Kelley [Ref1].

Picard iteration (2)

We cannot apply Picard iteration to nonlinear equations unless there is some special structure. For the commonly arising case $A(u)u=b(u)$ we can linearize the product $A(u)u$ to $A(u_{-})u$ and $b(u)$ as $b(u_{-})$. That is, we use the most previously computed approximation in $A$ and $b$ to arrive at a linear system for $u$:

$$A(u_{-})u = b(u_{-}){\thinspace .}$$

A relaxed iteration takes the form

$$A(u_{-})u^* = b(u_{-}),\quad u = \omega u^* + (1-\omega)u_{-}{\thinspace .}$$

In other words, we solve a system of nonlinear algebraic equations as a sequence of linear systems.

Algorithm for relaxed Picard iteration

Given $A(u)u=b(u)$ and an initial guess $u_{-}$, iterate until convergence:

1. solve $A(u_{-})u^* = b(u_{-})$ with respect to $u^*$
2. $u = \omega u^* + (1-\omega) u_{-}$
3. $u_{-}\ \leftarrow\ u$

Newton’s method (2)

The natural starting point for Newton’s method is the general nonlinear vector equation $F(u)=0$. As for a scalar equation, the idea is to approximate $F$ around a known value $u_{-}$ by a linear function $\hat F$, calculated from the first two terms of a Taylor expansion of $F$. In the multi-variate case these two terms become

$$F(u_{-}) + J(u_{-}) \cdot (u - u_{-}),$$

where $J$ is the Jacobian of $F$, defined by

$$J_{i,j} = \frac{\partial F_i}{\partial u_j}{\thinspace .}$$

So, the original nonlinear system is approximated by

$$\hat F(u) = F(u_{-}) + J(u_{-}) \cdot (u - u_{-})=0,$$

which is linear in $u$ and can be solved in a two-step procedure: first solve $J\delta u = -F(u_{-})$ with respect to the vector $\delta u$ and then update $u = u_{-} + \delta u$. A relaxation parameter can easily be incorporated:

$$u = \omega(u_{-} +\delta u) + (1-\omega)u_{-} = \omega_{-} + \omega\delta u{\thinspace .}$$

Algorithm for Newton’s method

Given $F(u)=0$ and an initial guess $u_{-}$, iterate until convergence:

1. solve $J\delta u = -F(u_{-})$ with respect to $\delta u$
2. $u = u_{-} + \omega)\delta u$
3. $u_{-}\ \leftarrow\ u$

For the special system with structure $A(u)u=b(u)$,

$$F_i = \sum_k A_{i,k}(u)u_k - b_i(u),$$

and

$$J_{i,j} = \sum_k \frac{\partial A_{i,k}}{\partial u_j}u_k + A_{i,j} - \frac{\partial b_i}{\partial u_j}{\thinspace .}$$

We realize that the Jacobian needed in Newton’s method consists of $A(u_{-})$ as in the Picard iteration plus two additional terms arising from the differentiation. Using the notation $A'(u)$ for $\partial A/\partial u$ (a quantity with three indices: $\partial A_{i,k}/\partial u_j$), and $b'(u)$ for $\partial b/\partial u$ (a quantity with two indices: $\partial b_i/\partial u_j$), we can write the linear system to be solved as

$$(A + A'u + b')\delta u = -Au + b,$$

or

$$(A(u_{-}) + A'(u_{-})u_{-} + b'(u_{i}))\delta u = -A(u_{-})u_{-} + b(u_{-}){\thinspace .}$$

Rearranging the terms demonstrates the difference from the system solved in each Picard iteration:

$$\underbrace{A(u_{-})(u_{-}+\delta u) - b(u_{-})}_{\hbox{Picard system}} +\, \gamma (A'(u_{-})u_{-} + b'(u_{i}))\delta u = 0{\thinspace .}$$

Here we have inserted a parameter $\gamma$ such that $\gamma=0$ gives the Picard system and $\gamma=1$ gives the Newton system. Such a parameter can be handy in software to easily switch between the methods.

Stopping criteria (2)

Let $||\cdot||$ be the standard Eucledian vector norm. Four termination criteria are much in use:

• Absolute change in solution: $||u - u_{-}||\leq \epsilon_u$
• Relative change in solution: $||u - u_{-}||\leq \epsilon_u ||u_0||$, where $u_0$ denotes the start value of $u_{-}$ in the iteration
• Absolute residual: $||F(u)|| \leq \epsilon_r$
• Relative residual: $||F(u)|| \leq \epsilon_r ||F(u_0)||$

To prevent divergent iterations to run forever, one terminates the iterations when the current number of iterations $k$ exceeds a maximum value $k_{\max}$.

The relative criteria are most used since they are not sensitive to the characteristic size of $u$. Nevertheless, the relative criteria can be misleading when the initial start value for the iteration is very close to the solution, since an unnecessary reduction in the error measure is enforced. In such cases the absolute criteria work better. It is common to combine the absolute and relative measures of the size of the residual, as in

$$||F(u)|| \leq \epsilon_{rr} ||F(u_0)|| + \epsilon_{ra},$$

where $\epsilon_{rr}$ is the tolerance in the relative criterion and $\epsilon_{ra}$ is the tolerance in the absolute criterion. With a very good initial guess for the iteration (typically the solution of a differential equation at the previous time level), the term $||F(u_0)||$ is small and $\epsilon_{ra}$ is the dominating tolerance. Otherwise, $\epsilon_{rr} ||F(u_0)||$ and the relative criterion dominates.

With the change in solution as criterion we can formulate and combined absolute and relative measure of the change in the solution:

$$||\delta u|| \leq \epsilon_{ur} ||u_0|| + \epsilon_{ua},$$

The ultimate termination criterion, combining the residual and the change in solution tests with a test on the maximum number of iterations allow, can be expressed as

$$\begin{split}||F(u)|| \leq \epsilon_{rr} ||F(u_0)|| + \epsilon_{ra} \hbox{ or } ||\delta u|| \leq \epsilon_{ur} ||u_0|| + \epsilon_{ua} \hbox{ or } k>k_{\max}{\thinspace .}\end{split}$$

Example: A nonlinear ODE model from epidemiology

The simplest model spreading of a disease, such as a flu, takes the form of a $2\times 2$ ODE system

$$S' = -\beta SI,$$

$$I' = \beta SI - \nu I,$$

where $S(t)$ is the number of people who can get ill (susceptibles) and $I(t)$ is the number of people who are ill (infected). The constants $\beta >0$ and $\nu >0$ must be given along with initial conditions $S(0)$ and $I(0)$.

Implicit time discretization

A Crank-Nicolson scheme leads to a $2\times 2$ system of nonlinear algebraic equations in the unknowns $S^{n+1}$ and $I^{n+1}$:

$$\frac{S^{n+1}-S^n}{\Delta t} = -\beta [SI]^{n+\frac{1}{2}} \approx -\frac{\beta}{2}(S^nI^n + S^{n+1}I^{n+1}),$$

$$\frac{I^{n+1}-I^n}{\Delta t} = \beta [SI]^{n+\frac{1}{2}} - \nu I^{n+\frac{1}{2}} \approx \frac{\beta}{2}(S^nI^n + S^{n+1}I^{n+1}) - \frac{\nu}{2}(I^n + I^{n+1}){\thinspace .}$$

Introducing $S$ for $S^{n+1}$, $S_1$ for $S^n$, $I$ for $I^{n+1}$, $I_1$ for $I^n$, we can rewrite the system as

(1) $$F_S(S,I) = S - S_1 + \frac{1}{2}\Delta t\beta(S_1I_1 + SI) = 0,$$

(2) $$F_I(S,I) = I - I_1 - \frac{1}{2}\Delta t\beta(S_1I_1 + SI) - \frac{1}{2}\Delta t\nu(I_1 + I) =0{\thinspace .}$$

A Picard iteration

We assume that we have approximations $S_{-}$ and $I_{-}$ to $S$ and $I$. A way of linearizing the only nonlinear term $SI$ is to write $I_{-}S$ in the $F_S=0$ equation and $S_{-}I$ in the $F_I=0$ equation, which also decouples the equations. Solving the resulting linear equations with respect to the unknowns $S$ and $I$ gives

$$\begin{split}S &= \frac{S_1 - \frac{1}{2}\Delta t\beta S_1I_1} {1 + \frac{1}{2}\Delta t\beta I_{-}}, \\ I &= \frac{I_1 + \frac{1}{2}\Delta t\beta S_1I_1} {1 - \frac{1}{2}\Delta t\beta S_{-} + \nu}{\thinspace .}\end{split}$$

The solutions $S$ and $I$ are stored in $S_{-}$ and $I_{-}$ and a new iteration is carried out.

Newton’s method (3)

The nonlinear system (1)-(2) can be written as $F(u)=0$ with $F=(F_S,F_I)$ and $u=(S,I)$. The Jacobian becomes

$$\begin{split}\renewcommand*{\arraystretch}{2} J = \left(\begin{array}{cc} \frac{\partial}{\partial S} F_S & \frac{\partial}{\partial I}F_S\\ \frac{\partial}{\partial S} F_I & \frac{\partial}{\partial I}F_I \end{array}\right) = \left(\begin{array}{cc} 1 + \frac{1}{2}\Delta t\beta I & \frac{1}{2}\Delta t\beta\\ - \frac{1}{2}\Delta t\beta S & 1 - \frac{1}{2}\Delta t\beta I - \frac{1}{2}\Delta t\nu \end{array}\right) {\thinspace .}\end{split}$$ $$The Newton system to be solved in each iteration is then$$ $$!bt$$ $$\begin{split}\renewcommand*{\arraystretch}{1.5} & \left(\begin{array}{cc} 1 + \frac{1}{2}\Delta t\beta I_{-} & \frac{1}{2}\Delta t\beta S_{-}\\ - \frac{1}{2}\Delta t\beta S_{-} & 1 - \frac{1}{2}\Delta t\beta I_{-} - \frac{1}{2}\Delta t\nu \end{array}\right) \left(\begin{array}{c} \delta S\\ \delta I \end{array}\right) =\\ & \qquad\qquad \left(\begin{array}{c} S_{-} - S_1 + \frac{1}{2}\Delta t\beta(S_1I_1 + S_{-}I_{-})\\ I_{-} - I_1 - \frac{1}{2}\Delta t\beta(S_1I_1 + S_{-}I_{-}) - \frac{1}{2}\Delta t\nu(I_1 + I_{-}) \end{array}\right)\end{split}$$

Remark. For this particular system explicit time integration methods work very well. The 4-th order Runge-Kutta method is an excellent balance between high accuracy, high efficiency, and simplicity.