Order and Error Constants

5.2. Order and Error Constants#

Theorem 5.2 (Global Error and Order of Accuracy)

If the local error associated to a multistep formula is \(\,O(h^{p+1})\) then, the global error is expected to be \(\,O(h^p)\), and we can say that the method is of order \(p\). To determine the order of a multistep method we can use the following approach.

The most general linear multistep method has the form

(5.1)#\[\begin{aligned} \sum_{i=0}^{k} \alpha_i y_{j+i} = h \sum_{i=0}^{k} \beta_i f_{j+i}, \hspace{1cm}\text{with}\ \alpha_k = 1 \end{aligned}\]

where \(k\) denotes the number of steps and \(\alpha_i\) and \(\beta_i\) are constants.

The exact solution of the differential equation at \(x = x + ih\), (i.e. \(\,y(x+ih)\,\)), also satisfies equation (5.1), so we can write

\[\begin{aligned} \sum_{i=0}^{k} \alpha_i ~ y{(x+ih)} ~=~ h \sum_{i=0}^{k} \beta_i ~ f{(x+ih)} ~+~ T \end{aligned}\]

Thus the local error \(T\) is given by

(5.2)#\[\begin{aligned} T ~=~ \sum_{i=0}^{k} \alpha_i ~ y{(x+ih)} ~-~ h \sum_{i=0}^{k} \beta_i ~ f{(x+ih)} \end{aligned}\]

Using the Taylor series expansion

\[\begin{aligned} y(x+ih) ~=~ y(x) + ihy'(x) + \frac{(ih)^2}{2!}y''(x) + \dots + \frac{(ih)^p}{p!}y^{(p)}(x)~, \end{aligned}\]

we can expand the right-hand side of equation (5.2) as:

  • \(\displaystyle i = 0\quad \alpha_0 y(x) - h\beta_0 y'(x) \)

  • \(\displaystyle i = 1\quad \alpha_1 \left[ y(x) + hy'(x) + \frac{h^2}{2!}y''(x) + \dots\right] -h\beta_1 \left[ y'(x) + hy''(x) + \frac{h^2}{2!}y'''(x) + \dots \right] \)

  • \(\displaystyle i = 2\quad \alpha_2 \left[ y(x) + 2hy'(x) + \frac{(2h)^2}{2!}y''(x) + \dots \right] -h\beta_2 \left[ y'(x) + 2hy''(x) + \frac{(2h)^2}{2!}y'''(x) + \dots \right]\)

  • \(\displaystyle i = 3\quad \alpha_3 \left[y(x) + 3hy'(x) + \frac{(3h)^2}{2!}y''(x) + \dots \right] -h\beta_3 \left[y'(x) + 3hy''(x) + \frac{(3h)^2}{2!}y'''(x) + \dots \right] \)

Thus, collecting the corresponding coefficients of \(y(x),\,y',\,...\) we have the following expressions [Note: we will add the stepsize \(h\), attached to each term, later in the equation for \(T\), shown below]

\[\begin{split}\begin{aligned} & \bigl[ (\alpha_0 + \alpha_1 + \alpha_2 + \dots + \alpha_k) \bigr]\,y(x) \\ & \bigl[ (\alpha_1 + 2\alpha_2 + 3\alpha_3 + \dots + k\alpha_k) - (\beta_0 + \beta_1 + \beta_2 + \dots + \beta_k) \bigr]\,y'(x) \\ & \bigl[ \frac{1}{2!}(\alpha_1 + 2{^2}\alpha_2 + 3{^2}\alpha_3 + \dots + k{^2}\alpha_k) - \frac{1}{1!}({\beta_1} + {2}\beta_2 + \dots + {k}\beta_k) \bigr]\,y''(x) \end{aligned}\end{split}\]

or in general,

\[\begin{aligned} \left[ \frac{1}{p!}(\alpha_1 + 2^p\alpha_2 + 3^p\alpha_3 + \dots + k^p\alpha_k) - \frac{1}{(p-1)!}(\beta_1 + 2^{p-1}\beta_2 + 3^{p-1}\beta_3 + \dots + k^{p-1}\beta_k) \right] y^{(p)}(x) \end{aligned}\]

for \(p = 2,3,\dots\)

Thus, let coefficients of \(y(x) = C_0\), \(~y'(x) = C_1\),…etc., then \(T\) can be written as,

\[\begin{aligned} T ~=~ C_0y(x) \,+\, C_1hy'(x) \,+\, C_2h^2y''(x) \,+\, \dots \,+\, C_p\,h^p\,y^{(p)}(x) \end{aligned}\]

The term \(C_{p+1}\,h^{p+1}\,y^{(p+1)}(x)~\) is called the principal or leading local truncation error and \(C_{p+1}\,\) is the error constant. For a polynomial of degree \(\,p\,\) we have \(~y^{(p+1)}(x) = 0~\) and hence the method is of order \(\,p\,\) if, \(~C_0 = C_1 = C_2 = \dots = C_p=0,~\) and \(~C_{p+1}\neq 0\).

Summary

  • \(C_0 = \alpha_0 + \alpha_1 + \alpha_2 + \cdots + \alpha_k\)

  • \(C_1 = (\alpha_1 + 2\alpha_2 + 3\alpha_3 + \dots + k\alpha_k) - (\beta_0 + \beta_1 + \beta_2 + \dots + \beta_k)\)

  • \(C_2 = \frac{1}{2!}(\alpha_1 + 2{^2}\alpha_2 + 3{^2}\alpha_3 + \dots + k{^2}\alpha_k) - \frac{1}{1!}({\beta_1} + {2}\beta_2 + \dots + {k}\beta_k)\)

  • For \(p=2, 3, \ldots\) :

    \[ C_p = \frac{1}{p!}(\alpha_1 + 2^p\alpha_2 + 3^p\alpha_3 + \dots + k^p\alpha_k) - \frac{1}{(p-1)!}(\beta_1 + 2^{p-1}\beta_2 + 3^{p-1}\beta_3 + \dots + k^{p-1}\beta_k) \]

Example 5.3

Find the order and the error constants of the following formulae,

  1. \(\displaystyle y_{j+4} = y_j + \frac{h}{3}(8f_{j+3} - 4f_{j+2} + 8f_{j+1})\)

    Solution

    A comparison with

    \[\begin{aligned} \sum_{i=0}^{4} \alpha_i y_{j+i} ~=~ h \sum_{i=0}^{4} \beta_i f_{j+i}, \end{aligned}\]

    gives,

    \[\begin{split}\begin{aligned} && \alpha_0 &= -1, & \alpha_1 &= 0, & \alpha_2 &= 0, & \alpha_3 &= 0, & \alpha_4 &= 1 && \\ && \beta_0 &= 0, & \beta_1 &= \frac{8}{3}, & \beta_2 &= -\frac{4}{3}, & \beta_3 &= \frac{8}{3}, & \beta_4 &= 0 && \end{aligned}\end{split}\]

    Thus

    \[\begin{split}\begin{aligned} C_0 ~&=~ \alpha_0 + \alpha_1 + \alpha_2 + \alpha_3 + \alpha_4 ~=~ 0 \phantom{\frac{}{0}}\\ C_1 ~&=~ (\alpha_1 + 2\alpha_2 + 3\alpha_3 + 4\alpha_4) - (\beta_0 + \beta_1 + \beta_2 + \beta_3 + \beta_4) ~=~ 0 \phantom{\frac{0}{0}}\\ C_2 ~&=~ \frac{1}{2!}(\alpha_1 + 2^2\alpha_2 + 3^2\alpha_3 + 4^2\alpha_4) - \frac{1}{1!}(\beta_1 + 2\beta_2 + 3\beta_3 + 4\beta_4) ~=~ 0\\ C_3 ~&=~ \frac{1}{3!}(\alpha_1 + 2^3\alpha_2 + 3^3\alpha_3 + 4^3\alpha_4) - \frac{1}{2!}(\beta_1 + 2^2\beta_2 + 3^2\beta_3 + 4^2\beta_4) ~=~ 0\\ C_4 ~&=~ \frac{1}{4!}(\alpha_1 + 2^4\alpha_2 + 3^4\alpha_3 + 4^4\alpha_4) - \frac{1}{3!}(\beta_1 + 2^3\beta_2 + 3^3\beta_3 + 4^3\beta_4) ~=~ 0\\ C_5 ~&=~ \frac{1}{5!}(\alpha_1 + 2^5\alpha_2 + 3^5\alpha_3 + 4^5\alpha_4) - \frac{1}{4!}(\beta_1 + 2^4\beta_2 + 3^4\beta_3 + 4^4\beta_4) ~=~ \frac{28}{90}=\frac{14}{45}\\ \end{aligned}\end{split}\]

    \(\therefore\) The method is of order \(4\) and the error constant is \(\dfrac{14}{45}\).

    Matlab code for error constant calculation:

    %matlab code: calculating the error constant for chapter 3 example 3(i)
%show results in fractions
format rational

a0=-1; a1=0;  a2=0;  a3=0; a4=1;
b0=0; b1=8/3; b2=-4/3; b3=8/3;  b4=0;

c0= a0+a1+a2+a3+a4;
c1= (a1+2*a2+3*a3+4*a4)-(b0+b1+b2+b3+b4);
c2= (a1+2^2*a2+3^2*a3+4^2*a4)/factorial(2)-(b1+2*b2+3*b3+4*b4)/factorial(1);
c3= (a1+2^3*a2+3^3*a3+4^3*a4)/factorial(3)-(b1+2^2*b2+3^2*b3+4^2*b4)/factorial(2);
c4= (a1+2^4*a2+3^4*a3+4^4*a4)/factorial(4)-(b1+2^3*b2+3^3*b3+4^3*b4)/factorial(3);
c5= (a1+2^5*a2+3^5*a3+4^5*a4)/factorial(5)-(b1+2^4*b2+3^4*b3+4^4*b4)/factorial(4);

disp([c0 c1 c2 c3 c4 c5])

  2. \(\displaystyle y_{j+3} ~=~ y_{j+2} + \frac{h}{12}(23f_{j+2} - 16f_{j+1} + 5f_j)\)

    Solution

    A comparison with

    \[\begin{aligned} \sum_{i=0}^{3} \alpha_i y_{j+i} ~=~ h \sum_{i=0}^{3} \beta_i f_{j+i} \end{aligned}\]

    gives,

    \[\begin{split}\begin{aligned} && \alpha_0 &= 0, & \alpha_1 &= 0, & \alpha_2 &= -1, & \alpha_3 &= 1 && \\ && \beta_0 &= \frac{5}{12}, & \beta_1 &= -\frac{16}{12}, & \beta_2 &= \frac{23}{12}, & \beta_3 &= 0 && \end{aligned}\end{split}\]

    Thus

    \[\begin{split}\begin{aligned} C_0 ~&=~ \alpha_0 + \alpha_1 + \alpha_2 + \alpha_3 ~=~ 0 \phantom{\frac{}{0}}\\ C_1 ~&=~ (\alpha_1 + 2\alpha_2 + 3\alpha_3) - (\beta_0 + \beta_1 + \beta_2 + \beta_3) ~=~ 0 \phantom{\frac{0}{0}}\\ C_2 ~&=~ \frac{1}{2!}(\alpha_1 + 2^2\alpha_2 + 3^2\alpha_3) - \frac{1}{1!}(\beta_1 + 2\beta_2 + 3\beta_3) ~=~ 0\\ C_3 ~&=~ \frac{1}{3!}(\alpha_1 + 2^3\alpha_2 + 3^3\alpha_3) - \frac{1}{2!}(\beta_1 + 2^2\beta_2 + 3^2\beta_3) ~=~ 0\\ C_4 ~&=~ \frac{1}{4!}(\alpha_1 + 2^4\alpha_2 + 3^4\alpha_3) - \frac{1}{3!}(\beta_1 + 2^3\beta_2 + 3^3\beta_3) ~=~ \frac{9}{24}\\ \end{aligned}\end{split}\]

    \(\therefore\) The method is of order \(3\) and the error constant is \(\dfrac{3}{8}\).

    Matlab code for error constant calculation:

    %matlab code: calculating the error constant for chapter 3 example 3(ii)
%show results in fractions
format rational

a0=0; a1=0;  a2=-1;  a3=1; 
b0=5/12; b1=-16/12; b2=23/12; b3=0;  

c0 = a0+a1+a2+a3;
c1 = (a1+2*a2+3*a3)-(b0+b1+b2+b3);
c2 = (a1+2^2*a2+3^2*a3)/factorial(2)-(b1+2*b2+3*b3)/factorial(1);
c3 = (a1+2^3*a2+3^3*a3)/factorial(3)-(b1+2^2*b2+3^2*b3)/factorial(2);
c4 = (a1+2^4*a2+3^4*a3)/factorial(4)-(b1+2^3*b2+3^3*b3)/factorial(3);
c5 = (a1+2^5*a2+3^5*a3)/factorial(5)-(b1+2^4*b2+3^4*b3)/factorial(4);

disp([c0 c1 c2 c3 c4 c5])

Note

The formulae for \(\,C_0\), \(C_1\), …, \(C_p\,\) can be used to derive a linear multistep method of given structure and maximal order (see Example 5.4, below).

Example 5.4

Construct an implicit linear two-step method of maximal order, containing one free parameter, and find its order.

Solution

For a two-step method \(\,k \!=\! 2;\ \alpha_2 \!=\! 1\,\), by hypothesis. Let \(\,\alpha_0 \!=\! a\,\) be the free parameter. There remain four undetermined coefficients \(\,\alpha_1,\, \beta_0,\, \beta_1\,\) and \(\,\beta_2\,\), and thus we can set \(\,C_0\!=\!C_1\!=\!C_2\!=\!C_3\!=\!0\,\).

Using the formulae for \(\,C_0,\,C_1,\,C_2\,\) and \(\,C_3\,\) we get

\[\begin{split}\begin{aligned} C_0 ~&=~ a + \alpha_1 + 1 ~=~ 0\\ C_1 ~&=~ \alpha_1 + 2 - (\beta_0 + \beta_1 + \beta_2) ~=~ 0\\ C_2 ~&=~ \frac{1}{2!}(\alpha_1 + 4) - (\beta_1 + 2\beta_2) ~=~ 0\\ C_3 ~&=~ \frac{1}{3!}(\alpha_1 + 8) - \frac{1}{2!}(\beta_1 + 4\beta_2) ~=~ 0 \end{aligned}\end{split}\]

Solving this set of equations gives

\[\begin{aligned} \alpha_1 &= -1 - a, & \beta_0 &= -\frac{1}{12}(1 + 5a), & \beta_1 &= \frac{2}{3}(1 - a), & \beta_2 &= \frac{1}{12}(5 + a) \end{aligned}\]

and the method is

(5.3)#\[\begin{aligned} y_{j+2} - (1 + a)y_{j+1} + ay_j ~=~ \frac{{h}}{12} \bigl[ (5 + a)f_{j+2} + 8(1 - a)f_{j+1} - (1 + 5a)f_j \bigr] \end{aligned}\]

Furthermore,

\[\begin{split}\begin{aligned} C_4 ~&=~ \frac{1}{4!}(\alpha_1 + 16) - \frac{1}{3!}(\beta_1 + 8\beta_2) ~=~ -\frac{1}{4!}(1 + a)\\ C_5 ~&=~ \frac{1}{5!}(\alpha_1 + 32) - \frac{1}{4!}(\beta_1 + 16\beta_2) ~=~ -\frac{1}{3\times 5!}(17 + 13a) \end{aligned}\end{split}\]

If \(\,a \!\neq\! -1\,\), then \(\,C_4 \!\neq\! 0\,\), and method (5.3) is of order \(3\). If \(\,a \!=\! -1\,\), then \(\,C_4 \!=\! 0\,,~ C_5 \!\neq\! 0\,\,\), and method (5.3), which is now Simpson’s rule, is of order \(4\).

Note that when \(\,a \!=\! 0\,\), then method (5.3) is the two-step Adams-Moulton method, and if \(\,a \!=\! -5\,\), it is an explicit method.

The following code has been tested under Matlab 2024a

syms a0 a1 b0 b1 b2

a2=1;

C0 = a0+a1+a2;
C1 = (a1+2*a2)-(b0+b1+b2);
C2 = (a1+2^2*a2)/factorial(2)-(b1+2*b2)/factorial(1);
C3 = (a1+2^3*a2)/factorial(3)-(b1+2^2*b2)/factorial(2);
C4 = (a1+2^4*a2)/factorial(4)-(b1+2^3*b2)/factorial(3);

solve([C0, C1, C2, C3, C4], [a0,a1,b0,b1,b2]) 

import math
import sympy as sp

a2=1;

a0, a1, b0, b1, b2 = sp.symbols('a0, a1, b0, b1, b2');

C0 = a0+a1+a2;
C1 = (a1+2*a2)-(b0+b1+b2);
C2 = (a1+2**2*a2)/math.factorial(2)-(b1+2*b2)/math.factorial(1);
C3 = (a1+2**3*a2)/math.factorial(3)-(b1+2**2*b2)/math.factorial(2);
C4 = (a1+2**4*a2)/math.factorial(4)-(b1+2**3*b2)/math.factorial(3);

solution=sp.solve([C0, C1, C2, C3, C4])

print(solution)

5.2.1. Computer programs for calculating error constants#

The following code has been tested under Matlab 2024a

 1%Matlab code for calculating the error constants of a linear multistep method
 2%
 3%  a0*y0+a1*y1+a2*y2+a3*y3+a4*y4=h*(b0*f0+b1*f1+b2*f2+b3*f3+b4*f4)
 4%
 5
 6format rational
 7
 8a0=-1; a1=0;  a2=0;  a3=0; a4=1;
 9b0=0; b1=8/3; b2=-4/3; b3=8/3;  b4=0;
10
11c0= a0+a1+a2+a3+a4;
12c1= (a1+2*a2+3*a3+4*a4)-(b0+b1+b2+b3+b4);
13c2= (a1+2^2*a2+3^2*a3+4^2*a4)/factorial(2)-(b1+2*b2+3*b3+4*b4)/factorial(1);
14c3= (a1+2^3*a2+3^3*a3+4^3*a4)/factorial(3)-(b1+2^2*b2+3^2*b3+4^2*b4)/factorial(2);
15c4= (a1+2^4*a2+3^4*a3+4^4*a4)/factorial(4)-(b1+2^3*b2+3^3*b3+4^3*b4)/factorial(3);
16c5= (a1+2^5*a2+3^5*a3+4^5*a4)/factorial(5)-(b1+2^4*b2+3^4*b3+4^4*b4)/factorial(4);
17
18disp([c0 c1 c2 c3 c4 c5])

 1"""
 2Python code for calculating the error constants of a linear multistep method
 3
 4a0*y0+a1*y1+a2*y2+a3*y3+a4*y4=h*(b0*f0+b1*f1+b2*f2+b3*f3+b4*f4)
 5"""
 6
 7
 8import math
 9from fractions import Fraction
10
11a0=Fraction( -1, 1);      
12a1=Fraction(  0, 1);      
13a2=Fraction(  0, 1);       
14a3=Fraction(  0, 1);     
15a4=Fraction(  1, 1);
16
17b0=Fraction(  0, 3);  
18b1=Fraction(  8, 3);  
19b2=Fraction( -4, 3);  
20b3=Fraction(  8, 3);  
21b4=Fraction(  0, 3);
22
23c0= a0+a1+a2+a3+a4;
24c1= (a1+2*a2+3*a3+4*a4)-(b0+b1+b2+b3+b4);
25c2= (a1+2**2*a2+3**2*a3+4**2*a4)/math.factorial(2)-(b1+2*b2+3*b3+4*b4)/math.factorial(1);
26c3= (a1+2**3*a2+3**3*a3+4**3*a4)/math.factorial(3)-(b1+2**2*b2+3**2*b3+4**2*b4)/math.factorial(2);
27c4= (a1+2**4*a2+3**4*a3+4**4*a4)/math.factorial(4)-(b1+2**3*b2+3**3*b3+4**3*b4)/math.factorial(3);
28c5= (a1+2**5*a2+3**5*a3+4**5*a4)/math.factorial(5)-(b1+2**4*b2+3**4*b3+4**4*b4)/math.factorial(4);
29
30print(c0, c1, c2, c3, c4, c5)
Contents