4.5.4 Various stop criteria
We let \( \delta y_i \) denote the increment between two consecutive iterations: $$ \begin{equation} \tag{4.115} \begin{array}{c} \delta y_i = y_i^{m+1} - y_i^m \ , \ i = 1,2,...,N, \qquad m = 0,1,...\\ \end{array} \end{equation} $$
where \( N \) denotes the number of physical grid points. Two classes of stop criteria based on absolute and relative metrics, may then be formulated as outlined in the following.
Stop criteria based on absolute values. $$ \begin{equation} \max(|\delta y_i|) < \varepsilon_a \tag{4.116} \end{equation} $$ $$ \begin{equation} \frac{1}{n}\sum_{i=1}^N |\delta y_i| < \varepsilon_a \tag{4.117} \end{equation} $$ $$ \begin{equation} \frac{1}{n}\sqrt{\sum_{i=1}^N(\delta y_i)^2} < \varepsilon_a \tag{4.118} \end{equation} $$
where \( \varepsilon_a \) is a suitable small number.
Stop criteria based on relative values $$ \begin{equation} \max \left(\left| \frac{\delta y_i}{y_i^{m+1}}\right| \right) < \varepsilon_r , \qquad y_i^{m+1} \not= 0 \tag{4.119} \end{equation} $$ $$ \begin{equation} \frac{\sum_{i=1}^N|\delta y_i|}{\sum_{i=1}^N\left|y_i^{m+1}\right|} < \varepsilon_r \tag{4.120} \end{equation} $$ $$ \begin{equation} \frac{\max(|\delta y_i|)}{\max(|y_i^{m+1}|)} < \varepsilon_r \tag{4.121} \end{equation} $$
Note that if the quantity of interest, i.e. \( y_i \), is of order 1, the absolute and the relative criteria are equivalent. Normally we will prefer a relative criterion as these comply with the expression number of correct digits. If cases where the values of \( y_i \) are small in the whole computational domain (i.e. for all i) we may resort to a absolute value criterion.
4.5.5 Example: Usage of the various stop criteria
In this example we will demonstrate the use of the various stop criteria outlined in the previous section on the problem in section (3.2 Shooting methods for boundary value problems with nonlinear ODEs), where we will compute the second solution \( y_{II}) \) as illustrated in Figure (37).
For convenience we reformulate the ODE (4.92) : $$ \begin{equation} y''(x) = \frac{3}{2}y^2 \tag{4.122} \end{equation} $$
with the boundary conditions: $$ \begin{equation} y(0) = 4 , \qquad y(1) = 1 \tag{4.123} \end{equation} $$
We let our initial guess be expressed by the parabola \( y_s = 20(x-x^2) \) , \( 0 < x < 1 \), which is plotted along with the solution \( y_{II} \) , in Figure (61).
Figure 61: Initial parabolic guess and solution \( y_{II} \).
For this example it is natural to use a relative stop criterion as we know that the solution has values in the range \( [4, -10.68] \). After nine iterations the following table results:
| Iter.nr. | \( t_{r2} \) | \( t_{r3} \) |
| 1 | \( 7.25\cdot 10^{-1} \) | \( 7.55\cdot 10^{-1} \) |
| 2 | \( 8.70\cdot 10^{-1} \) | \( 8.20\cdot 10^{-1} \) |
| 3 | \( 8.49\cdot 10^{-1} \) | \( 6.65\cdot 10^{-1} \) |
| 4 | \( 4.96\cdot 10^{-1} \) | \( 5.85\cdot 10^{-1} \) |
| 5 | \( 2.10\cdot 10^{-1} \) | \( 2.80\cdot 10^{-1} \) |
| 6 | \( 4.62\cdot 10^{-2} \) | \( 6.07\cdot 10^{-2} \) |
| 7 | \( 2.24\cdot 10^{-3} \) | \( 2.70\cdot 10^{-3} \) |
| 8 | \( 4.68\cdot 10^{-6} \) | \( 5.79\cdot 10^{-6} \) |
| 9 | \( 2.26\cdot 10^{-11} \) | \( 2.55\cdot 10^{-11} \) |
The iteration evolution is representative for Newton-algorithms as the initial values are far away from the correct solution. The errors decreases slowly for the first six iterations, but afterwards we observe quadratic convergence. The two criteria do not differ significantly for this case.
A python code using a relative criteria is shown in avvikr below:
# src-ch3/avvikr.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)
# numerical parameters
h = 0.05 # steplength.
n = round(1./h)-1 # number of equations. h should be set so that n is an integer
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
itmax = 15
for it in range(itmax):
"""iteration process of linearized system of equations"""
# 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.ones(n)*fac*0.5*ym**2)
d[n-1] = -1.
d[0] = -4.
ym1 = tdma(a,b,c,d) # solution
dy = np.abs(ym1-ym)
tr2 = np.max(dy)/np.max(np.abs(ym1))
tr3 = np.sum(dy)/np.sum(np.abs(ym1))
ym = ym1
print('it = {}, tr2 = {}, tr3 = {} '.format(it, tr2, tr3))
ya = 4./(1+x)**2
feil = np.abs((ym1-ya)/ya)
# plotting:
legends=[] # empty list to append legends as plots are generated
# Add the labels
plot(x,ym1,"b") # plot numerical solution
legends.append('y2 numerical')
plot(x,ya,"r") # plot analytical solution
legends.append('y1 analytic')
legend(legends,loc='best',frameon=False) # Add the legends
ylabel('y')
xlabel('x')
grid(b=True, which='both', color='0.65',linestyle='-')
show()
A code where the absolute criteria is used is almost identical, except for the substitution of \( t_{r2} \) and \( t_{r3} \)
ta1 = max(dy);
ta2 = sum(dy)/n;
ta3 = sqrt(dot(dy,dy))/n;
# src-ch3/avvika.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)
# numerical parameters
h = 0.05 # steplength.
n = round(1./h)-1 # number of equations. h should be set so that n is an integer
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
itmax = 15
for it in range(itmax):
"""iteration process of linearized system of equations"""
# 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.ones(n)*fac*0.5*ym**2)
d[n-1] = -1.
d[0] = -4.
ym1 = tdma(a,b,c,d) # solution
dy = np.abs(ym1-ym)
ta1 = np.max(dy) # abs stop critertia 1
ta2 = np.sum(dy)/n # abs stop critertia 2
ta3 = np.sqrt(np.dot(dy,dy))/n # abs stop critertia 3
ym = ym1
print('it = {}, ta2 = {}, ta3 = {} '.format(it, ta2, ta3))
xv = np.linspace(h,1-h,n)
ya = 4./(1+xv)**2
# plotting:
legends=[] # empty list to append legends as plots are generated
# Add the labels
plot(xv,ym1,"b")
legends.append('y2 numerical')
#plot(xv,ya,"r")
#legends.append('y1 analytic')
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()