$$ \newcommand{\D}{\displaystyle} \renewcommand{\eqref}[1]{Eq.~(\ref{#1})} $$
Numerical Methods for Engineers

 

 

 

4.5 Linearization of nonlinear algebraic equations

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)

Figure 60: Discretization one-dimensional boundary value problem on a unit interval.

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.

4.5.1 Picard 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)))

The method of Picard linearization. The advantage with the method of Picard linearization, is that the linearization process is very simple. However, the simplicity comes at a prize, since it converges slowly and often requires a relatively good initial guess.

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.

4.5.2 Newton linearization

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.

4.5.3 Example: Newton linearization of various nonlinear ODEs

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.