In this section we will present methods for linearization of nonlinear algebraic equations and will use the problem in section (4.4.2 Cooling rib with variable cross-section) as an illustrative example: $$ \begin{equation} y''(x) = \frac{3}{2}y^2 \tag{4.92} \end{equation} $$

with the following boundary conditions: $$ \begin{equation} y(0) = 4, \qquad y(1) = 1 \tag{4.93} \end{equation} $$

From (4.4.2 Cooling rib with variable cross-section) we know that one of the two analytical solutions is given by: $$ \begin{equation} y = \frac{4}{(1+x)^2} \tag{4.94} \end{equation} $$

We use central differences for the term \( y''(x) \) to discretize equation (4.92): $$ \begin{equation*} \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2}= \frac{3}{2} \; y_i^2 \end{equation*} $$

which after collection of terms may be presented as: $$ \begin{equation} y_{i-1} - \left( 2+\frac{3}{2}h^2 \; y_i \right) \; y_i + y_{i+1} = 0 \tag{4.95} \end{equation} $$

We have discretized the unit interval \( [0, 1] \) in \( N \) parts with a corresponding mesh size \( h=1/N \) and discrete coordinates denoted by \( x_i = h\cdot i \) for \( i = 0,1,\ldots, N+1 \) (see Fig 60)

By substituting the known boundary values \( y_0 = 4 \) and \( y_{N+1} = 1 \) in equation (4.95) we get the following system of equations: $$ \begin{align} -\left( 2+\frac{3}{2}h^2 \; y_1 \right) y_1 + y_2 &= -4 \nonumber \\ & \vdots \nonumber \\ y_{i-1} - \left(2+\frac{3}{2}h^2 \; y_i\right) y_i+y_{i+1} &= 0 \tag{4.96}\\ & \vdots \nonumber\\ y_{N-1} - \left( 2 + \frac{3}{2}h^2 \; y_N\right) y_N &= -1 \nonumber \end{align} $$

where \( i = 1,2,\ldots, N-1 \). The coefficient matrix is tridiagonal, but the system is nonlinear. And as we have no direct analytical solutions for such nonlinear algebraic systems, we must linearize the system and develop an iteration strategy which hopefully will converge to a solution to the nonlinear problem. In the following we will present two methods for linearization.

Due to the nonlinearities in equation (4.95), we must develop an iteration strategy and to do so we let \( y_i^{m+1} \) and \( y_i^m \) be the solution of the discretized equation (4.96) at iterations \( m+1 \) and \( m \), respectively.

With Picard linearization we linearize the nonlinear terms simply by replacement of dependent terms at iteration \( m+1 \) with corresponding terms at iteration \( m \) until only linear terms are left. In equation (4.96) the nonlinearity is caused by the \( y^2 \)-term and the corresponding terms need to be linearized.

By using the iteration nomenclature, equation (4.96) takes the form: $$ \begin{align} -\left( 2 + \frac{3}{2}y_1^{m+1}h^2\right) \; y_1^{m+1} + y_2^{m+1} &= -4 \nonumber \\ & \vdots \nonumber \\ y_{i-1}^{m+1} - \left( 2 + \frac{3}{2}y_i^{m+1}h^2\right) \; y_i^{m+1} + y_{i+1}^{m+1} &= 0 \tag{4.97}\\ & \vdots \nonumber \\ y_{N-1}^{m+1} - \left( 2 + \frac{3}{2}y_N^{m+1}h^2\right) \; y_N^{m+1} &= -1 \nonumber \end{align} $$

where \( i = 2,3,...N-1 \), and the iteration counter \( m =
0,1,2,\ldots \). The linearization is performed by replacing
\( \frac{3}{2}h^2y_1^{m+1} \) , \( \frac{3}{2}h^2y_i^{m+1} \) and
\( \frac{3}{2}h^2y_N^{m+1} \) by \( \frac{3}{2}h^2y_1^m \) ,
\( \frac{3}{2}h^2y_i^m \) and \( \frac{3}{2}h^2y_N^m \) such that the
following *linear system* is obtained:
$$
\begin{align}
-\left( 2 + \frac{3}{2}y_1^{m}h^2\right) \; y_1^{m+1} + y_2^{m+1} &= -4 \nonumber \\
& \vdots \nonumber\\
y_{i-1}^{m+1} - \left( 2 + \frac{3}{2}y_i^{m}h^2\right)\; y_i^{m+1} + y_{i+1}^{m+1} &= 0 \tag{4.98} \\
& \vdots \nonumber\\
y_{N-1}^{m+1} - \left( 2 + \frac{3}{2}y_N^{m}h^2\right) \; y_N^{m+1} &= -1 \nonumber
\end{align}
$$

We have with the procedure above discretized our analytical problem in
equation (4.95) to a *linear*, tridiagonal system which may
readily be solved with sparse solvers like the ones we introduced in 4.4.1 Example: Heat exchanger with constant cross-section.

Note that all values at iteration level *m*, i.e. \( y_i^{m} \) are
*known* and the only unknowns in (4.98) are \( y_i^{m-1} \). Thus,
to start the iteration procedure we need an initial value for
\( y_i^0 \). Unless a well informed initial guess is available, it is
common to start with \( y_i^0 = 0 \), \( i = 1,2,...N \). Naturally, initial
values *close* to the numerical solution will result in faster
convergence.

Additionally, we must test for *diagonal dominance* of the iteration
scheme, which for equation (4.98) means:
$$
\begin{equation}
\left| 2+\frac{3}{2}y_i^mh^2 \right| \ge 2 \ , \ i=1,2,...,N
\tag{4.99}
\end{equation}
$$

We observe from equation (4.99) that we have a diagonally
dominant iteration matrix when \( y_i^m > 0 \). By inspection the
analytical solutions in Figure 37, we observe that this
condition is fulfilled for only one of the two solutions. In cases
when the condition of diagonal dominance is *not* fulfilled, make sure
your solver allows for pivotation, e.g. like `scipy.sparse.linalg`

as
demonstrated in 4.4.1 Example: Heat exchanger with constant cross-section. However, for demonstration purposes we
have in this example implemented a version of the classical
tri-diagonal solver `tdma`

along with a version with pivotation
`tripiv`

in the module `TRIdiagonalSolvers`

.

In cases where the one has no experience with the iteration process, e.g. while implementing, one may find it beneficial to use the number of iterations as a stopping criteria, rather than an error/residual based stop criterion. However, after you have gained some experience with the particular problem at hand you may use other stop criteria as demonstrated in the section 4.5.4 Various stop criteria.

In the python-code `delay34.py`

below we have implemented the
procedure of Picard linearization as outlined above.

```
# src-ch3/delay34.py;TRIdiagonalSolvers.py @ git@lrhgit/tkt4140/src/src-ch3/TRIdiagonalSolvers.py;
import numpy as np
from matplotlib.pyplot import *
from TRIdiagonalSolvers import tdma
from TRIdiagonalSolvers import tripiv
# change some default values to make plots more readable
LNWDT=3; FNT=11
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT
font = {'size' : 16}; rc('font', **font)
h = 0.05 # steplength
n = int(round(1./h)-1) # number of equations
fac = (3./2.)*h**2
Nmax = 30 # max number of iterations
# iterative solution with Nmax iterations
#initilizing:
ym = np.zeros(n) # initial guess of y. ym = np.zeros(n) will converge to y1. ym = -20*x*(1-x) for instance will converge to y2
for m in range(Nmax):
"""iteration process of linearized system of equations using delay method"""
# need to set a, b, c and d vector inside for loop since indices are changed in tdma solver
a = np.ones(n)
c = np.ones(n)
b = -(np.ones(n)*2. + fac*ym)
d = np.zeros(n)
d[n-1] = -1.
d[0] = -4.
ym1 = tdma(a,b,c,d) # solution
if(m>0):
max_dym = np.max(np.abs((ym1 - ym)/ym))
print(('m = {} \t max dym = {:6.4g} '.format(m, max_dym)))
ym = ym1
# compute analytical solution
xv = np.linspace(h, 1 - h, n)
ya = 4./(1 + xv)**2
error = np.abs((ym1 - ya)/ya)
print(('\nThe max relative error after {} iterations is: {:4.3g}'.format(Nmax,np.max(error))))
```

The maximal relative deviation between consecutive iterations is
stored in the variable `max_dym`

. We observe that `max_dym`

decreases
for each iteration as the solution slowly converges to the analytical
solution \( y_I = \dfrac{4}{(1+x)^2} \) which was demonstrated in (37) in section (3.2 Shooting methods for boundary value problems with nonlinear ODEs).

After `Nmax`

iterations we compute the `error`

by :

```
error = np.abs((ym1 - ya)/ya)
```

Note that `error`

is a vector by construction, since `ym1`

and `ya`

are vectors. Thus, to print a scalar measure of the error after \tt{Nmax}
iterations we issue the following command:

```
print('\nThe max relative error after {} iterations is: {:4.3g}'.format(Nmax,np.max(error)))
```

Once we have gained some experience on how the iteration procedure develops, we may make use of a stop criterion based on relative changes in the iterative solutions as shown in the python-snippet below:

```
# iterative solution with max_dym stop criterion
m=0; RelTol=1.0e-8
max_dym=1.0
ym = np.zeros(n) # initial guess of y. ym = np.zeros(n) will converge to y1. ym = -20*x*(1-x) for instance will converge to y2
a = np.ones(n)
c = np.ones(n)
while (m<Nmax and max_dym>RelTol):
d = np.zeros(n)
d[n-1] = -1.
d[0] = -4.
b = -(np.ones(n)*2. + fac*ym)
ym1 = tdma(a,b,c,d) # solution
if(m>0):
max_dym = np.max(np.abs((ym1 - ym)/ym))
print(('m = {} \t max dym = {:6.4g} '.format(m, max_dym)))
m+=1
ym = ym1
error = np.abs((ym1 - ya)/ya)
print(('\nThe max relative after {} iterations is: {:4.3g}'.format(m,np.max(error))))
# print results nicely with pandas
print_results=0
if (print_results):
import pandas as pd
data=np.column_stack((xv,ym1,ya))
data_labeled = pd.DataFrame(data, columns=['xv','ym','ya'])
print((data_labeled.round(3).to_string(index=False)))
# plotting:
legends=[] # empty list to append legends as plots are generated
# Add the labels
plot(xv,ya,"r") # plot analytical solution
legends.append('analytic')
plot(xv,ym1,"b--") # plot numerical solution
legends.append('delayed linearization')
legend(legends,loc='best',frameon=False) # Add the legends
ylabel('y')
xlabel('x')
#grid(b=True, which='both', axis='both',linestyle='-')
grid(b=True, which='both', color='0.65',linestyle='-')
show()
```

Notice that we must set `max_dym>RelTol`

prior to the `while-loop`

to start the iterative process.

You may experiment with the code above to find that for a fixed mesh
size `h=0.05`

, the relative error will not be smaller than \( \approx
5.6 \cdot 10^{-4} \), regardless of how low you set `RelTol`

, i.e. after
a converged solution is obtained for the given mesh size. Naturally,
the error will be reduced if you reduce the mesh size `h`

. Observe
also that the number of iterations needed to obtain a converged
solution for a given mesh size of `h=0.05`

is smaller than the
`Nmax=30`

, which was our original upper limit.

**Motivating example of linearization**

Let us look at a motivating example which is simple to use for nonlinearities in products. For convenience we introduce \( \delta y_i \) as the increment between two consecutive \( y_i \)-values: $$ \begin{equation} y^{m+1}_i = y_i^m + \delta y_i \tag{4.100} \end{equation} $$

By assuming \( \delta y_i \) to be small, the nonlinear term in equation (4.95) may then be linearized using (4.100) in the following way: $$ \begin{equation} \left( y_i^{m+1} \right)^2 = \left( y_i^m + \delta y_i \right)^2 = (y_i^m)^2 + 2y_i^m \, \delta y_i + (\delta y_i^m)^2 \approx (y_i^m)^2 + 2 \, y_i^m \, \delta y_i = y_i^m \; \left (2 y_i^{m+1}-y_i^m \right ) \tag{4.101} \end{equation} $$

That is, quadratic expression in equation (4.101) is linearized by neglecting the quadratic terms \( (\delta y)^2 \) which are small compared to the other terms in (4.101).

Substitution of the linearization (4.101) in (4.96) results in the following system of equations: $$ \begin{equation} \begin{array}{rcl} -(2+3y_1^m h^2) \; y_1^{m+1} + y_2^{m+1} &= -\frac{3}{2} \, (y_i^mh)^2 -4\\ \vdots\\ y_{i-1}^{m+1} - (2+3y_i^mh^2) \; y_i^{m+1} + y_{i+1}^{m+1}&= -\frac{3}{2} \, (y_i^mh)^2\\ \vdots\\ y_{N-1}^{m+1}- (2+3y_N^mh^2) \; y_N^{m+1} &= -\frac{3}{2} \, (y_N^mh)^2 -1\\ \end{array} \tag{4.102} \end{equation} $$

where \( i = 1,2,\ldots,N-1 \) and \( m=0,1,\ldots \).

The resulting equation system (4.102), is a tridiagonal system which may be solved with the tdma-solver as for the previous example. An implementation is given below:

```
# src-ch3/taylor34.py;TRIdiagonalSolvers.py @ git@lrhgit/tkt4140/src/src-ch3/TRIdiagonalSolvers.py;
import numpy as np
from matplotlib.pyplot import *
from TRIdiagonalSolvers import tdma
from TRIdiagonalSolvers import tripiv
# change some default values to make plots more readable
LNWDT=3; FNT=11
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT
font = {'size' : 16}; rc('font', **font)
h = 0.05 # steplength
n = int(round(1./h)-1) # number of equations
fac = (3.)*h**2
#initilizing:
x = np.linspace(h,1-h,n)
ym = -20*x*(1-x) # initial guess of y.
#ym = np.zeros(n) # will converge to y1. ym = -20*x*(1-x) for instance will converge to y2
it, itmax, dymax, RelTol = 0, 15, 1., 10**-10
legends=[]
while (dymax > RelTol) and (it < itmax):
"""iteration process of linearized system of equations using taylor
"""
plot(x,ym) # plot ym for iteration No. it
legends.append(str(it))
it = it + 1
a = np.ones(n)
c = np.ones(n)
b = -(np.ones(n)*2. + fac*ym)
d = -(fac*0.5*ym**2)
d[n-1] = d[n-1]-1
d[0] = d[0]-4
ym1 = tdma(a,b,c,d) # solution
dymax = np.max(np.abs((ym1-ym)/ym))
ym = ym1
print(('it = {}, dymax = {} '.format(it, dymax)))
legend(legends,loc='best',frameon=False)
ya = 4./(1+x)**2
feil = np.abs((ym1-ya)/ya)
print("\n")
for l in range(len(x)):
print('x = {}, y = {}, ya = {}'.format(x[l], ym1[l], ya[l]))
figure()
# plotting:
legends=[] # empty list to append legends as plots are generated
# Add the labels
plot(x,ya,"r") # plot analytical solution
legends.append('y1 analytic')
plot(x,ym1,"b--") # plot numerical solution
legends.append('y2 taylor linearization')
legend(legends,loc='best',frameon=False) # Add the legends
ylabel('y')
xlabel('x')
#grid(b=True, which='both', axis='both',linestyle='-')
grid(b=True, which='both', color='0.65',linestyle='-')
show()
```

You may experiment with by changing the initial, guessed solution, and observe what happens with the resulting converged solution. Which solution does the numerical solution converge to?

Observe that the iteration process converges much faster than for the Picard linearization approach. Notably, somewhat slow in the first iterations, but after some iterations quadratic convergence is obtained.

**Generalization: Newton linearization**

The linearization in (4.100) and (4.101) may seen as
truncated Taylor expansions of the dependent variable at iteration
*m*, where only the first two terms are retained.

To generalize, we let *F* denote the nonlinear term to be linearized
and assume *F* to be a function of the *dependent* variable *z* at the
discrete point \( x_i \) (or *i* for short) at iteration *m*.

We will linearize \( F(z_i)\big|^{m+1}\equiv F(z_i)_{m+1} \), and use the
latter notation in the following. A Taylor expansion may then be
expressed at the point *i* and iteration *m*
$$
\begin{equation}
F(z_i)_{m+1} \approx F(z_i)_m + \delta z_i \left(\dfrac{\partial F}{dz_i}\right)_m
\tag{4.103}
\end{equation}
$$

where \( \delta z_i = z_i^{m+1} - z_i^m \).

Our simple example in equation (4.101) may be rendered in this slightly more generic representation: $$ \begin{equation*} \delta z_i \to \delta y_i, \qquad z_i^{m+1} \to y_i^{m+1}, \qquad z_i^m \to y_i^m, \qquad F(y_i)_{m+1} = (y_i^{m+1})^2 \ , \ F(y_i)_m = (y_i^m)^2 \end{equation*} $$

which after substitution in (4.103) yields: \( (y_i^{m+1})^2 \approx (y_i^m)^2 + 2y_i^m \cdot \delta y_i \) which is identical to what was achieved in equation (4.101) .

The linearization in equation (4.103) is often referred to as
*Newton-linearization*.

In many applications, the nonlinear term is composed \( z \) at various discrete, spatial locations and this has to be accounted for in the linearization. Assume to begin with that we have terms involving both \( z_i \) and \( z_{i+1} \). The Taylor expansion will then take the form: $$ \begin{equation} F(z_i,z_{i+1})_{m+1} \approx F(z_i,z_{i+1})_m + \left(\dfrac{\partial F}{\partial z_i}\right)_m \; \delta z_i + \left( \dfrac{\partial F}{\partial z_{i+1}} \right)_m \; \delta z_{i+1} \tag{4.104} \end{equation} $$

where: $$ \begin{equation*} \delta z_i = z_i^{m+1} - z_i^m , \qquad \ \delta z_{i+1} = z_{i+1}^{m+1} - z_{i+1}^m \end{equation*} $$

When three spatial indices, e.g. \( z_{i-1},z_i \) and \( z_{i+1} \), are involved, the Taylor expansion becomes: $$ \begin{equation} F(z_{i-1},z_i,z_{i+1})_{m+1} \approx F(z_{i-1},z_i,z_{i+1})_m + \left( \dfrac{\partial F}{\partial z_{i-1}}\right)_m \; \delta z_{i-1} + \left( \dfrac{\partial F}{\partial z_{i}}\right)_m \; \delta z_{i} + \left( \dfrac{\partial F}{\partial z_{i+1}}\right)_m \; \delta z_{i+1} \tag{4.105} \end{equation} $$

where $$ \begin{equation*} \delta z_{i-1} = z_{i-1}^{m+1} - z_{i-1}^m , \qquad \ \delta z_{i} = z_{i}^{m+1} - z_{i}^m , \qquad \ \delta z_{i+1} = z_{i+1}^{m+1} - z_{i+1}^m \label{} \end{equation*} $$

Naturally, the same procedure may be expanded whenever the nonlinear term is composed of instances of the dependent variable at more spatial locations.

**Nonlinear ODE 1**

In this example of a nonlinear ODE we consider: $$ \begin{equation} y''(x) + y(x)\sqrt{y(x)} = 0 \label{} \end{equation} $$

which after discretization with central differences may be presented: $$ \begin{equation} y_{i+1}^{m+1} - 2y_i^{m+1} + y_{i-1}^{m+1} + h^2y_i^{m+1}\sqrt{y_i^{m+1}}= 0 \tag{4.106} \end{equation} $$

In equation (4.106) it is the nonlinear term \( y_i^{m+1} \; \sqrt{y_i^{m+1}} \) which must be linearized. As we in this example only have one index we make use of the simplest Taylor expansion (4.103) with \( z_i \to y_i \) and \( F(y_i) = y_i^{\frac{3}{2}} \) : $$ \begin{equation} F(y_i)_{m+1} \approx F(y_i)_m + \delta y_i \left( \dfrac{\partial F}{\partial y_i}\right)_m = (y_i^m)^{\frac{3}{2}} + \delta y_i \; \frac{3}{2} \; (y_i^m)^{\frac{1}{2}} \tag{4.107} \end{equation} $$ where \( \delta y_i = y_i^{m+1} - y_i^m \). Equivalently by using the squareroot-symbols: $$ \begin{equation} y_i^{m+1}\sqrt{y_i^{m+1}}\approx y_i^m\sqrt{y_i^m}+\frac{3}{2}\; \sqrt{y_i^m} \; \delta y_i \tag{4.108} \end{equation} $$

Observe that we have succeeded in the linearization process as all occurrences of \( y_i^{m+1} \) in the two equivalent discretization in equation (4.107) and (4.108) are of power one.

After substitution in equation (4.106) the following *linear* difference equation results:
$$
\begin{equation} \tag{4.109}
y_{i-1}^{m+1} +\left( \frac{3}{2} h^2 \sqrt{y_i^m}-2\right) y_i^{m+1} + y_{i+1}^{m+1} = \frac{h^2}{2}y_i^m\sqrt{y_i^m}
\end{equation}
$$

**Nonlinear ODE 2**

As a second example of a nonlinear ODE consider: $$ \begin{equation} y''(x) + \sin\left (y(x)\right) = 0 \label{} \end{equation} $$

which after discretization by central differences may be rendered: $$ \begin{equation} y_{i+1}^{m+1} - 2y_i^{m+1} + y_{i-1}^{m+1} + h^2 \sin \left( y_i^{m+1} \right) = 0 \tag{4.110} \end{equation} $$

In equation (4.110) we identify \( \sin \left( y_i^{m+1} \right) \) as the nonlinear term in need of linearization. Again, we have only one index and may make use of (4.103) with \( z_i \to y_i \) and \( F(y_i) = \sin (y_i) \): $$ \begin{equation*} F(y_i)_{m+1} = \sin(y_i^{m+1}) \approx F(y_i)_m + \delta y_i \left( \dfrac{\partial F}{\partial y_i}\right)_m = \sin(y_i^m) + \delta y_i \; \cos(y_i^m) \end{equation*} $$

which after substitution in equation (4.110) leads to the following difference equation: $$ \begin{equation} \tag{4.111} y_{i-1}^{m+1} + (h^2 \cos(y_i^m) - 2) y_i^{m+1} + y_{i+1}^{m+1} = h^2(y_i^m \cos(y_i^m) - \sin(y_i^m)) \end{equation} $$

**Nonlinear ODE 3**

Consider the ODE \( y''(x) + y(x)\sqrt{y'(x)}= 0 \) which after discretization with central differences may be presented: $$ \begin{equation} \tag{4.112} y_{i+1}^{m+1} -2y_i^{m+1} + y_{i-1}^{m+1} + \alpha \; y_i^{m+1}\sqrt{y_{i+1}^{m+1} - y_{i-1}^{m+1}} = 0 \ ,\text{ der $\alpha = h\sqrt{h/2}$} \end{equation} $$

The term \( y_i^{m+1}\sqrt{y_{i+1}^{m+1} - y_{i-1}^{m+1}} \) is nonlinear and must be linearized. As three indices are involved in the term we utilize (4.105) with \( z_{i-1} \to y_{i-1} \), \( \ z_i \to y_i \), \( z_{i+1} \to y_{i+1} \), and finally \( F(y_{i-1},y_i,y_{i+1})_m = y_i \sqrt{y_{i+1}-y_{i-1}} \).

We present the individual terms in equation(4.105) : $$ \begin{equation*} \begin{array}{ll} F(y_{i-1},y_i,y_{i+1})_m = y_i^m\sqrt{y_{i+1}^m- y_{i-1}^m}, & \D \left( \frac{\partial F}{\partial y_{i-1}}\right)_m = - \frac{y_i^m}{2\sqrt{y_{i+1}^m - y_{i-1}^m}}\\ \D \left( \frac{\partial F}{\partial y_i}\right)_m = \sqrt{y_{i+1}^m - y_{i-1}^m}, & \D \left( \frac{\partial F}{\partial y_{i+1}}\right)_m = \frac{y_i^m}{2\sqrt{y_{i+1}^m - y_{i-1}^m}} \end{array} \end{equation*} $$

By collecting terms we get: $$ \begin{equation} y_i^{m+1} \sqrt{y_{i+1}^{m+1} - y_{i-1}^{m+1}} \approx y_i^m \sqrt{y_{i+1}^m - y_{i-1}^m} - \frac{y_i^m}{2 \sqrt{y_{i+1}^m - y_{i-1}^m}} \; \delta y_{i-1} + \sqrt{y_{i+1}^m - y_{i-1}^m} \; \delta y_i + \frac{y_i^m}{2 \sqrt{y_{i+1}^m - y_{i-1}^m}} \; \delta y_{i+1} \tag{4.113} \end{equation} $$

Note that equation (4.113) is linear as \( y_{i-1}^{m+1} \), \( y_i^{m+1} \) and \( y_{i+1}^{m+1} \) are only in first power. By substitution in equation (4.112) we get the linear discretized equation: $$ \begin{equation} \tag{4.114} \begin{array}{c} \left( 1 - \alpha\frac{y_i^m}{2g^m}\right) \; y_{i-1}^{m+1} - (2 - \alpha g^m)\; y_i^{m+1} + \left( 1 + \alpha\frac{y_i^m}{2g^m}\right) \; y_{i+1}^{m+1}\\\\ = \alpha\frac{y_i^m}{2g^m}(y_{i+1}^m-y_{i-1}^m)\\\\ \text{der } g^m = \sqrt{y_{i+1}^m-y_{i-1}^m} \end{array} \end{equation} $$

In the next section we will show how we can operate directly with the incremental quantities \( \delta y_{i-1}^{m+1} \), \( \delta y_{i}^{m+1} \) and \( \delta y_{i+1}^{m+1} \).

**Comparison with the method of Picard linearization**

How would the examples above be presented with the method of Picard linearization?

Equation (4.106) takes the following form the Picard linearization: $$ \begin{equation} y_{i-1}^{m+1}+(h^2\sqrt{y_i^{m+1}}-2)\; y_i^{m+1}+y_{i-1}^{m+1} = 0 \label{} \end{equation} $$

The coefficient for the \( y_i^{m+1} \)-term contains \( \sqrt{y_i^{m+1}} \) which will be replaced by \( \sqrt{y_i^m} \) such that: $$ \begin{equation} y_{i-1}^{m+1}+(h^2\sqrt{y_i^m}-2)\; y_i^{m+1}+y_{i-1}^{m+1} = 0 \label{} \end{equation} $$

For Picard linearization equation (4.110) become: $$ \begin{equation} y_{i-1}^{m+1} - 2y_i^{m+1}+y_{i+1}^{m+1} = -h^2\sin(y_i^m)=0 \label{} \end{equation} $$

whereas equation (4.112) looks like: $$ \begin{equation} y_{i-1}^{m+1} + (\alpha \; \sqrt{y_{i+1}^{m+1}-y_{i-1}^{m+1}} - 2) \; y_i^{m+1}+y_{i+1}^{m+1} = 0 \label{} \end{equation} $$

with Picard linearization. The coefficient in front of the \( y_i^{m+1} \)-term contains \( \sqrt{y_{i+1}^{m+1}-y_{i-1}^{m+1}} \) which is replaced by \( \sqrt{y_{i+1}^m-y_{i-1}^m} \) such that the equation becomes: $$ \begin{equation*} y_{i-1}^{m+1} + (\alpha\; \sqrt{y_{i+1}^{m} - y_{i-1}^{m}}-2)\; y_i^{m+1} + y_{i+1}^{m+1} = 0 \end{equation*} $$

When we use the method of Picard linearization, the system first has to be presented in a form which identifies the coefficients in the equation system. If the coefficients have dependent variable at iteration \( m+1 \), the are to be replaced with the corresponding value at iteration \( m \). Often this procedure is identical to a Taylor expansion truncated after the first term.

Thus, we may expect that the method of Picard linearization will converge more slowly compared to methods which use more terms in the Taylor expansion. The method of Picard linearization is most frequently used for partial differential equations where e.g. the material parameters are functions of the dependent variables, such as temperature dependent heat conduction.