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: $$ \begin{equation} y_i^{m+1} = y_i^m + \delta y_i , \qquad i=1,2,\ldots,N \tag{4.127} \end{equation} $$
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()