
## 4.5.6 Linearizations of nonlinear equations on incremental form

As seen previously in equation (4.101), we may introduce the increment $$\delta y_i = y_i^{m+1} - y_i^m$$ in a nonlinear term as exemplified by: \begin{align} (y_i^{m+1})^2 &= (y_i^m + \delta y_i)^2 = (y_i^m)^2 + 2y_i^m \; \delta y_i + (\delta y_i)^2 \tag{4.124}\\ &\approx (y_i^m)^2 + 2y_i^m \; \delta y_i = y_i^m(2y_i^{m+1} - y_i^m) \tag{4.125} \end{align}

i.e. we linearize by ignoring higher order terms of the increments, and solve equation (4.125) with respect to $$y_i^{m+1}$$.

Alternatively, we may think of $$\delta y_i$$ as the unknown quantity and solve the resulting equation system with respect to the increments. For example we may substitute $$y_k^{m+1} = y_k^m + \delta y_k$$ for $$k = i-1, \; i, \; i+1$$ and $$(y_i^{m+1})^2 \approx (y_i^m)^2 + 2 \, y_i^m \; \delta y_i$$ into equation (4.97) (or equivalently (4.102)), the following equation system results: \begin{align} -(2+3y_1^mh^2)\;\delta y_1 + \delta y_2 &= -(y_2^m - 2y_1^m + 4) + \frac{3}{2}\; (y_1^mh)^2 \nonumber \\ \delta y_{i-1} - (2+3y_i^mh^2)\; \delta y_i + \delta y_{i+1} &= -(y_{i+1}^m - 2y_i^m + y_{i-1}^m) + \frac{3}{2} \; (y_i^mh)^2 \tag{4.126}\\ \delta y_{N-1} - (2+3y_N^mh^2)\;\delta y_N &= -(1-2y_N^m + y_{N-1}^m)+\frac{3}{2}\;(y_N^mh)^2 \nonumber \end{align}

where $$i=2,3,\ldots,N-1$$ and $$m=0,1,2,\ldots$$. In the equation system above we have incorporated the boundary conditions $$y_0 = 4$$, $$\delta y_0 = 0$$, $$y_{N+1} = 1$$, and $$\delta y_{N+1}=0$$. For each iteration we update the y-values by: $$$$y_i^{m+1} = y_i^m + \delta y_i , \qquad i=1,2,\ldots,N \tag{4.127}$$$$

After having gained some experience with how the iteration proceed, we may use stop criteria such ad $$\max|\delta y_i| < \varepsilon_1$$ or $$\max|\delta y_i / y_i^{m+1}| < \varepsilon_2$$, $$i=1,2,\ldots,N$$.

We have implemented the incremental formulation in the code below.

# src-ch3/delta34.py;ODEschemes.py @ git@lrhgit/tkt4140/src/src-ch3/ODEschemes.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
nitr = 6 # number of 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

it, itmax, dymax, relTol = 0, 10, 1., 10**-5 # numerical tollerance limits

while (dymax > relTol) and (it < itmax):
"""iteration process of linearized system of equations"""
it = it + 1
a = np.ones(n)
c = np.ones(n)
b = -(np.ones(n)*2. + fac*ym)
d = (np.ones(n)*fac*0.5*ym**2)
for j in range(1,n-1):
d[j] = d[j] - (ym[j+1] - 2*ym[j] +ym[j-1])

d[n-1] = d[n-1]-(1 - 2*ym[n-1] + ym[n-2])
d[0] = d[0] -(ym[1] - 2*ym[0] + 4)

dy = tdma(a,b,c,d) # solution
ym = ym + dy
dymax  = np.max(np.abs((dy)/ym))

print('it = {},  dymax = {} '.format(it, dymax))

xv = np.linspace(h,1-h,n)
ya = 4./(1+xv)**2

print("\n")

for l in range(len(xv)):
print('x = {},  y = {},  ya = {}'.format(xv[l], ym[l], ya[l]))

legends=[] # empty list to append legends as plots are generated

plot(xv,ya,"r") # plot analytical solution
legends.append('y1 analytic')
plot(xv,ym,"b--") # plot numerical solution
legends.append('delta linearization')
# Add the labels
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()