Preface

This section provides basic estimates for used numerical methods.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
Return to Part III of the course APMA0330

Glossary

Error estimates

Now that we have learned how to use the numerical methods to find a numerical solution to simple first-order differential equations, we need to develop some practical guidelines to help us estimate the accuracy of various methods. Because we have replaced a differential equation by a difference equation, our numerical solution cannot be identically equal, in general, to the “true” solution of the original differential equation. In general, the discrepancy between the two solutions has two causes. One cause is that computers do not store numbers with infinite precision, but rather to a maximum number of digits that is hardware and software dependent. Most programming languages allow the programmer to distinguish between floating point or real numbers, that is, numbers with decimal points, and integer or fixed point numbers. Arithmetic with numbers represented by integers is exact. However in general, we cannot solve a differential equation using integer arithmetic. Hence, arithmetic operations such as addition and division, which involve real numbers, can introduce an error, called the roundoff error. For example, if a computer only stored real numbers to two significant figures, the product 2.1×3.2 would be stored as 6.7 rather than the exact value 6.72. The ignificance of roundoff errors is that they accumulate as the number of mathematical operations increases. Ideally we should choose algorithms that do not significantly magnify the roundoff error, for example, we should avoid subtracting numbers that are nearly the same in magnitude.

The other source of the discrepancy between the true answer and the computed answer is the error associated with the choice of algorithm. This error is called the truncation error. A truncation error would exist even on an idealized computer that stored floating point numbers with infinite precision and hence had no roundoff error. Because this error depends on the choice of algorithm and hence can be controlled by the programmer, you should be motivated to learn more about numerical analysis and the estimation of truncation errors. However, there is no general prescription for the best method for obtaining numerical solutions of differential equations. We will find in later sections that each method has advantages and disadvantages and the proper selection depends on the nature of the solution, which you might not know in advance, and on your objectives. How accurate must the answer be? Over how large an interval do you need the solution? What kind of computer are you using? How much computer time and personal time do you have?

In practice, we usually can determine the accuracy of our numerical solution by reducing the value of ∆t until the numerical solution is unchanged at the desired level of accuracy. Of course, we have to be careful not to choose ∆t too small, because too many steps would be required and the total computation time and roundoff error would increase.

In addition to accuracy, another important consideration is the stability of an algorithm. For example, it might happen that the numerical results are very good for short times, but diverge from the true solution for longer times. This divergence might occur if small errors in the algorithm are multiplied many times, causing the error to grow geometrically. Such an algorithm is said to be unstable for the particular problem. One way to determine the accuracy of a numerical solution is to repeat the calculation with a smaller step size and compare the results. If the two calculations agree to p decimal places, we can reasonably assume that the results are correct to p decimal places.

We present first the Mean Value Theorem from calculus that have an impact on numerical methods. A special case of this theorem was first described by Parameshvara (1370--1460), from the Kerala school of astronomy and mathematics in India. A restricted form of the theorem was proved by the French mathematician Michel Rolle (1652--1719) in 1691; the result was what is now known as Rolle's theorem, and was proved only for polynomials, without the techniques of calculus. The mean value theorem in its modern form was stated and proved by the French mathematician Baron Augustin-Louis Cauchy (1789--1857) in 1823.

Theorem (Mean Value Theorem): Let \( f:\, [a,b] \,\mapsto \,\mathbb{R} \) be a continuous function on the closed interval [𝑎,,b], and differentiable on the open interval (𝑎,b), where 𝑎<b. Then there exists some ξ in (𝑎,b) such that

\[ f(b) - f(a) = f' (\xi ) \left( b-a \right) . \qquad\qquad⧫ \]

Theorem (Cauchy Mean Value Theorem): If functions f and g are both continuous on the closed interval [a,b], and differentiable on the open interval (𝑎,b), then there exists some \( \xi \in (a,b) , \) such that

\[ \left( f(b) - f(a) \right) g' (\xi ) = f' (\xi ) \left( g(b) - g(a) \right) . \qquad\qquad⧫ \]

Theorem (Mean Value Theorem for definite integrals): If \( f:\, [a,b] \,\mapsto \,\mathbb{R} \) is integrable and \( g:\, [a,b] \,\mapsto \,\mathbb{R} \) is an integrable function that does not change sign on [a, b], then there exists some \( \xi \in (a,b) , \) such that

\[ \int_a^b g(t)\,f(t)\,{\text d} t = f(\xi ) \int_a^b g(t) \,{\text d} t . \qquad\qquad⧫ \]

 

Error Estimates

We discussed so far different numerical algorithms to solve the initial value problems for first order differential equations:

\[ y' = f\left( x, y \right) , \qquad y(x_0 )= y_0 , \]
where the slope function f is given and defined on some compact 2D domain to which the initial point \( \left( x_0 , y_0 \right) \) belongs. We associate with this problem a uniform grid (only for simplicity, in practice it is almost always variable) of points \( x_n = x_0 + n\,h , \quad n=0,1,2,,\ldots , \) where h is the step size. We always assume that the given initial value problem has a unique solution, which we denote by \( y= \phi (x) . \) A numerical approximation leads to the sequence of ordinates \( y_n , \quad n=0,1,2,\ldots , \) that we wish approximate the true values \( \phi (x_n ) . \) Their difference \( E_n = \phi (x_n ) -y_n \) is known as the global truncation error.

The methods we introduce for approximating the solution of an initial value problem are called difference methods or discrete variable methods. The solution is approximated at a set of discrete points called a grid (or mesh) of points. An elementary single-step method has the form \( y_{n=1} = y_n + h\, \Phi (x_n, y_n ) \) for some function Φ called an increment function.

The errors incurred in a single integration step are of two types: 1. The local truncation error or discretization error - the error introduced by the approximation of the differential equation by a difference equation. 2. Errors due to the deviation of the numerical solution from the exact theoretical solution of the difference equation. Included in this class are round-off errors, due to the inability of evaluating real numbers with infinite precision with the use of computer (i.e., computers usually operate on fixed word length) , and errors which are incurred if the difference equation is implicit and is not solved exactly at each step.
If a multistep method is used, an additional source of error results from the use of an auxiliary method (usually a single-step method, e.g., Runge--Kutta method), to develop the needed starting values for the multistep method.

In the numerical solution of an ODE, a sequence of approximations yn to the actual solution \( y = \phi (x ) \) is generated. Roughly speaking, the stability of a numerical method refers to the behavior of the difference or error, \( y_n - \phi (x_n ) , \) as n becomes large. To make our exposition simple, we consider the propagation of error in a one-step method by studying the particular representative one-step method of Euler:

\[ y_{n+1} = y_{n} + h\, f\left( x_n, y_n \right) , \quad n=0,1,2,\ldots . \]
This can be done by determining the relation of the error at step n+1 to the error at step n. To do this let \( \phi_n = \phi (x_n ) \) denote the true value at x = xn of the actual solution \( y = \phi (x) \) to the initial value problem \( y' = f(x,y) , \quad y(x_0 )= y_0 . \) For simplicity, we consider only the uniform grid of points \( x_n = x_0 + n\,h , \quad n=0,1,2,\ldots , \) with a fixed step size h. Then the total accumulated solution error εn at step n is defined by
\[ \varepsilon_n = \phi_n - y_n , \qquad n=0,1,2,\ldots . \]
The numerical values computed by Euler's algorithm satisfy the relation
\[ y_{n+1} = y_{n} + h\, f\left( x_n, y_n \right) - R_{n+1} , \]
where Rn+1, denotes the round-off errors resulting from evaluating Euler's rule. Similarly the true values of actual solution satisfy the relation
\[ \phi_{n+1} = \phi_{n} + h\, f\left( x_n, \phi_n \right) + T_{n+1} , \]
where Tn+1 denotes the local truncation error.

Subtracting the latter from the former yields

\[ \phi_{n+1} - y_{n+1} = \phi_{n} - y_n + h\left[ f\left( x_n, \phi_n \right) - f\left( x_n, y_n \right) \right] + E_{n+1} , \]
where \( E_{n+1} = T_{n+1} + R_{n+1} . \) This yields
\[ \varepsilon_{n+1} = \varepsilon_{n} + h\left[ f\left( x_n, \phi_n \right) - f\left( x_n, y_n \right) \right] + E_{n+1} . \]
Applying the mean value theorem to the terms inside the bracket, this relation between successive errors can be written as
\[ \varepsilon_{n+1} = \varepsilon_{n} + h\left( \phi_n - y_n \right) f_y \left( x_n, \eta_n \right) + E_{n+1} , \]
where fy denotes \( \frac{\partial f}{\partial y} \) and ηn lies between yn and \( \phi_n . \) If \( |f_y (x , y )| \le K \) and \( |E_{n+1} | \le E , \) where K and E are both positive constants, the error expression above can be replaced by a related first-order difference equation
\[ |\varepsilon_{n+1} | \le | \varepsilon_{n} | \left[ 1 + h\,K \right] + E . \]
Now if \( |\varepsilon_{n} | \le | e_n | , \) it follows that
\[ |\varepsilon_{n+1} | \le | e_{n} | \left[ 1 + h\,K \right] + E . \]
The propagated error is bounded by the solution of the related first-order difference equation
\[ e_{n+1} = e_{n} \, q + E ,\qquad \mbox{with the initial condition} \quad e_0 =0 , \]
where \( q = 1 + h\,K . \) This means that the error propagates linearly. If \( |q| = |1 + h\,K| > 1 , \) the error increases while for \( |q| = |1 + h\,K| < 1 , \) decreases. This observation leads to the condition:
\[ -2 < h\, K < 0 \qquad \mbox{or} \qquad \partial f/\partial y <0 . \]
The solution of this first-order difference equation for en can be found by successive substitution as follows since it is actually a geometric series:
\[ e_{n} = \frac{1- q^n}{1-q} \, E . \]
It follows then that the propagated error in Euler's algorithm is bounded by the expression
\[ |\varepsilon_{n+1} | \le \frac{\left( 1 + h\,K \right)^{n+1} -1}{h\,K} \, E . \]
Now we can make the following conclusions.

 

The global discretization error ek is defined by
\[ e_k = \phi (x_k ) - y_k , \qquad k=0,1,2,\ldots , m . \]
It is the difference between the unique solution and the solution obtained by the discrete variable method.

The local discretization error εk+1 is defined by

\[ \varepsilon_{k+1} = \phi (x_{k+1} ) - y_k - h\, \Phi (x_k , y_k ), \qquad k=0,1,2,\ldots , m-1 . \]
It is the error committed in the single step from xk to xk+1.

Theorem (Precision of Euler's Method): Assume that \( y= \phi (x) \) is the solution to the initial value problem

\[ y' = f(x,y) \qquad\mbox{over interval} \quad [a,b]=[x_0 , x_m ] \quad\mbox{subject to} \quad y(x_0 )= y_0 . \]
If \( \phi (x) \) has two continuous derivatives on the interval [𝑎,b] and \( \left\{ \left( x_k , y_k \right) \right\}_{k=0}^m \) is the sequence of approximations generated by Euler's method, then
\[ \begin{split} |e_k | &= \left\vert \phi( x_k ) - y_k \right\vert = O(h) , \\ | \varepsilon_{k+1} | &= \left\vert \phi( x_{k+1} ) - y_k - h\, f(x_k , y_k ) \right\vert = O\left( h^2 \right) . \end{split} \]
The error at the end b = xm of the interval [𝑎,b] is called the final global error:
\[ E(\phi (b) , h) = \left\vert \phi (b) - y_m \right\vert = O(h) . \qquad\qquad⧫ \]
Suppose that the slope function f(x,y) is continuously differentiable with respect to both variables, x and y. Suppose also that the solution y = ϕ(x) of the initial value problem y' = f(x,y), y(x0) = y0 exist and unique on teh interval 𝑎 ≤ x ≤ b. WE are going to show that therte exists a positive constant M so that the global truncation error for Euler's rule satisfies
\[ E(\phi (b) , h) = \left\vert \phi (b) - y_m \right\vert \le M\,h . \]
Thus, Euler's rule is a first-order numerical method. Since f(x,y) is continuously differentiable, the solution y(x) has two continuous derivatives that is equal to

\[ y'' = f_x (x,y) + y'\,f_y (x,y) = f_x (x,y) + f(x,y)\,f_y (x,y) . \]
Now we estimate the local error using xn+1 = xn + h:
\begin{align*} e_{n+1} &= y\left( x_{n+1} \right) - y_{n+1} = y\left( x_{n} + h \right) - y_{n+1} \\ &= y \left( x_{n} \right) + h\,y' \left( x_{n} \right) + R_1 (x_n , h) - \left[ y_n + h\,f(x_n , y_n ) \right] \\ &= y \left( x_{n} \right) - y_n + h \left[ y' \left( x_{n} \right) - f \left( x_n , y_n \right) \right] + R_1 (x_n , h) \\ &= e_n + h \left[ y' \left( x_{n} \right) - f \left( x_n , y_n \right) \right] + R_1 (x_n , h) \end{align*}
Now we estimate the second and the third term in the above equation. Using the fact that y(x) is a solution of the differential equation y' = f(x,y), the second term is estimated as
\[ \left\vert y' \left( x_{n} \right) - f \left( x_n , y_n \right) \right\vert = \left\vert f \left( x_n , y (x_n ) \right) - f \left( x_n , y_n \right) \right\vert . \]
Applying the mean-value theorem, we get
\[ \left\vert y' \left( x_{n} \right) - f \left( x_n , y_n \right) \right\vert \le M_1 \left\vert e_n \right\vert , \]
where M1 is an upper bound for \( \left\vert f_y (x,y) \right\vert , \) taken over both x and y. Since the remainder in Taylor's approximation is
\[ R_1 (x,h) = \frac{1}{2}\, y'' (\xi )\, h^2 , \]
for some (unknown) point ξ, we have
\[ \left\vert R_1 (x,h) \right\vert \le M_2 h^2, \]
where 2 M2 is the maximum of \( \left\vert y'' (x,y) \right\vert , \) on the interval [𝑎,b]. Now taking the absolute value of the local error leads to
\[ \left\vert e_{n+1} \right\vert \le \left\vert e_{n} \right\vert + h\,M_1 \left\vert e_{n} \right\vert + h^2 M_2 . \]
Let &eps;n be the solution of the difference equation
\[ \varepsilon_{n+1} = \varepsilon_n + h\,M_1 \left\vert e_{n} \right\vert + h^2 M_2 , \qquad \varepsilon_0 = \left\vert e_{0} \right\vert . \]
It is not hard to show by induction that
\[ \varepsilon_n \ge \left\vert e_n \right\vert \qquad\mbox{for all $n$}. \]
Hence, to estimate \( \left\vert e_n \right\vert , \) it is sufficies to solve the first order difference equation \( \varepsilon_{n+1} = \varepsilon_n + h\,M_1 \left\vert e_{n} \right\vert + h^2 M_2 : \)
\[ \varepsilon_n = \left( 1 + h\,M_1 \right)^n \varepsilon_0 + \left[ \sum_{i=0}^{n-1} \left( 1 + h\,M_1 \right)^i \right] h^2 M_2 , \]
which is verified by substitution. Since &epsilonn ≤ &epsilonn+1, it follows that the largest &epsilonn is the last one &epsilonN. This yields
\[ E(h) \le \varepsilon_N . \]
Therefore, we need to estimate
\[ \varepsilon_N = \left( 1 + h\,M_1 \right)^N \varepsilon_0 + \left[ \sum_{i=0}^{N-1} \left( 1 + h\,M_1 \right)^i \right] h^2 M_2 . \]
The first term is
\[ \left( 1 + h\,M_1 \right)^N = \left( 1 + h\,M_1 \right)^{(b-a)/h} , \]
which converges to \( e^{M_1 (b-a)} \) because of
\[ \lim_{h \to 0} \left( 1 + h\,\alpha \right)^{\beta /h} = e^{\alpha \beta} . \]
The second term can be estimated using the formula of geometric progression
\[ \sum_{i=0}^{n-1} u^i = \frac{1-u^n}{1-u} . \]
Then
\begin{align*} \left[ \sum_{i=0}^{N-1} \left( 1 + h\,M_1 \right)^i \right] h^2 M_2 &= \left[ \frac{1- (1+ h\,M_1 )^N}{1 - (1+h\,M_1 )}\right] h^2 M_2 \\ &= \frac{(1+h\,M_1 )^{(b-a)/h} -1}{h\,M_1} \,h^2 M_2 \\ & \le \left[ e^{M_1 (b-a)} -1 \right] h \left( \frac{M_2}{M_1} \right) \\ & \le e^{M_1 (b-a)} h \left( \frac{M_2}{M_1} \right) . \end{align*}
This yields the final estimation
\[ \varepsilon_N \le e^{M_1 (b-a)} e_0 + e^{M_1 (b-a)} \left( \frac{M_2}{M_1} \right) h . \]
If e0 = 0, then
\[ \varepsilon_N \le M\,h , \]
where \( M = M_2 e^{M_1 (b-a)} / M_1 . \)    ▣
From the above theorem, we see that halving the step size has the effect of essentially halving the error. However, even then the accuracy was not as good as we probably would have liked. Of course we could just keep decreasing the step size (provided we did not take h to be so small that round-off errors started to play a role) to increase the accuracy, but then the number of steps we would have to take would make the calculations very cumbersome. A better approach is to derive methods that have a higher order of accuracy.

Example: Consider the initial value problem for the linear differential equation

\[ y' = x^2 y - 1.5 \,y , \qquad y(0) =1. \]
Its solution was found previously
DSolve[{y'[x] == x^2 *y[x] - 1.5*y[x], y[0] == 1}, y[x], x]
\[ y = \phi (x) = e^{x^3 /3 - 1.5\,x} . \]

We combine all steps into subroutine:
EulerMethod[n_,x0_,y0_,xf_,f_]:=Module[{Y,h,i},
h=(xf-x0)/n;
Y=y0;
For[i=0,i<n,i++,
Y=Y+f[x0+i*h,Y]*h]; Y]
In the above code,
n is number of steps
x0 = initial condition for x
y0 = initial condition for y
xf = final value of x at which y is desired.
f = slope function for the differential equation dy/dx = f(x,y)
nth=7;
Do[
Nn[i]=2^i;
H[i]=(xf-x0)/Nn[i];
AV[i]=EulerMethod[2^i,x0,y0,xf,f];
Et[i]=EV-AV[i];
et[i]=Abs[(Et[i]/EV)]*100.0;
If[i>0,
Ea[i]=AV[i]-AV[i-1];
ea[i]=Abs[Ea[i]/AV[i]]*100.0;
sig[i]=Floor[(2-Log[10,ea[i]/0.5])];
If[sig[i]<0,sig[i]=0];
] ,{i,0,nth}];
This loop calculates the following
AV = approximate value of the ODE solution using Euler's Method by calling the module EulerMethod
Et = true error
et = absolute relative true error percentage
Ea = approximate error
ea = absolute relative approximate error percentage
sig = least number of significant digits correct in approximation

The following code uses Euler's Method to calculate intermediate step values for the purpose of displaying the solution over the range specified. The number of steps used is the maximum value, n (in this example, n=128).

n=128;
X[0]=x0;
Y[0]=y0;
h=(xf-x0)/n;
For[i=1,i<=n,i++,
X[i]=x0+i*h;
Y[i]=Y[i-1]+f[X[i-1],Y[i-1]]*h; ];
data = Table[{X[i], Y[i]}, {i, 0, n}];
plot2 = ListPlot[data, Joined -> True,
PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 1]}];
Finally, we plot the exact solution along with Euler's approximation:
Show[plot, plot2, PlotLabel -> "Exact and Approximate Solution"]
Explicit solution (red) and its Euler's approximation (blue).

Now we compare the value at the final point (denoted by xf) with a step size

data = Table[{H[i], AV[i]}, {i, 0, nth}];
plot3 = ListPlot[data, Joined -> True,
PlotStyle -> {Thickness[0.006], RGBColor[0.5, 0.5, 0]}, PlotLabel ->
"Approximate value of the solution of the ODE\nat x = xf as a \ function of step size"]
Dependence of the value at the final point with step size.
data = Table[{Nn[i], AV[i]}, {i, 0, nth}];
plot4 = ListPlot[data, Joined -> True,
PlotStyle -> {Thickness[0.006], RGBColor[0, 0.5, 0.5]}, PlotLabel -> "Approximate value of the solution of the ODE \n at x = xf as a \ function of number of steps "]
Dependence of the value at the final point with the number of steps.
We plot dependences of errors with the step size.
soln = DSolve[{f[x, y[x]] == y'[x], y[x0] == y0}, y, x]
EV = (y /. First[soln])[xf]
data = Table[{Nn[i], Et[i]}, {i, 0, nth}]; plot5 = ListPlot[data, Joined -> True,
PlotStyle -> {Thickness[0.006], RGBColor[0.5, 0, 0.5]}, PlotLabel ->
"True error as a function of number of steps"]
True error as a function of number of steps.
data = Table[{Nn[i], et[i]}, {i, 0, nth}];
plot6 = ListPlot[data, AxesOrigin -> {0, 0},
Joined -> True, PlotStyle -> {Thickness[0.006], RGBColor[0.3, 0.3, 0.4]},
PlotLabel -> "Absolute relative true error percentage \n as a function of number of steps"]
Absolute relative true error percentage.
data = Table[{Nn[i], Ea[i]}, {i, 1, nth}];
plot7 = ListPlot[data, AxesOrigin -> {0, 0},
Joined -> True, PlotStyle -> {Thickness[0.006], RGBColor[0.7, 0.3, 0]},
PlotLabel -> "Approximate error \nas a function of number of steps"]
Approximate error.
data = Table[{Nn[i], ea[i]}, {i, 1, nth}];
plot9 = ListPlot[data, AxesOrigin -> {0, 0},
Joined -> True, PlotStyle -> {Thickness[0.006], RGBColor[0, 0.3, 0.7]},
PlotLabel -> "Absolute relative approximate error percentage \nas a function of \ number of steps"]
Absolute relative approximate error.
   ■

 

Theorem (Precision of Heun's Method): Assume that \( y= \phi (x) \) is the solution to the initial value problem

\[ y' = f(x,y) \qquad\mbox{over interval} \quad [a,b]=[x_0 , x_m ] \quad\mbox{subject to} \quad y(x_0 )= y_0 . \]
If \( \phi (x) \) has three continuous derivatives on the interval [𝑎,b] and \( \left\{ \left( x_k , y_k \right) \right\}_{k=0}^m \) is the sequence of approximations generated by Heun's method, then
\[ \begin{split} |e_k | &= \left\vert \phi( x_k ) - y_k \right\vert = O\left( h^2 \right) , \\ | \varepsilon_{k+1} | &= \left\vert \phi( x_{k+1} ) - y_k - h\, \Phi (x_k , y_k ) \right\vert = O\left( h^3 \right) , \end{split} \]
where \( \Phi (x_k , y_k ) = y_k + \left( \frac{h}{2} \right) \left[ f(x_k , y_k ) + f \left( x_{k+1} , y_k + h\,f(x_k ,y_k ) \right) \right] .\) The final global error at the end of the interval satisfies:
\[ E(\phi (b) , h) = \left\vert \phi (b) - y_m \right\vert = O\left( h^2 \right) . \qquad \qquad⧫ \]

 

Accuracy

Sometimes you may be interested to find out what methods are being used in NDSolve.
Here you can see the coefficients of the default 2(1) embedded pair.

NDSolve`EmbeddedExplicitRungeKuttaCoefficients[2, Infinity]
Out[9]=
{{{1}, {1/2, 1/2}}, {1/2, 1/2, 0}, {1, 1}, {-(1/2), 2/3, -(1/6)}}

You also may want to compare some of the different methods to see how they perform for a specific problem.

Implicit Runge--Kutta methods have a number of desirable properties.

The Gauss--Legendre methods, for example, are self-adjoint, meaning that they provide the same solution when integrating forward or backward in time.


http://reference.wolfram.com/language/tutorial/NDSolveImplicitRungeKutta.html

A generic framework for implicit Runge--Kutta methods has been implemented. The focus so far is on methods with interesting geometric properties and currently covers the following schemes:

"ImplicitRungeKuttaGaussCoefficients"
"ImplicitRungeKuttaLobattoIIIACoefficients"
"ImplicitRungeKuttaLobattoIIIBCoefficients"
"ImplicitRungeKuttaLobattoIIICCoefficients"
"ImplicitRungeKuttaRadauIACoefficients"
"ImplicitRungeKuttaRadauIIACoefficients"

The derivation of the method coefficients can be carried out to arbitrary order and arbitrary precision.

 

Errors in Numerical Approximations

The main question of any approximation is whether the numerical approximations \( y_1 , y_2 , y_3 , \ldots \) approach the corresponding values of the actual solution? We want to use a step size that is small enough to ensure the required accuracy. An unnecessarily small step length slows down the calculations and in some cases may even cause a loss of accuracy.

There are three fundamental sources of error in approximating the solution of an initial value problem numerically.

  1. The formula, or algorithm, used in the calculations is an approximating one.
  2. Except for the first step, the input data used in the calculations are only approximations to the actual values of the solution at the specified points.
  3. The computer used for the calculations has finite precision; in other words, at each stage only a finite number of digits can be retained.

Let us temporarily assume that our computer can execute all computations exactly; that is, it can retain infinitely many digits (if necessary) at each step. Then the difference En between the actual solution \( y = \phi (x) \) of the initial value problem \( y' = f(x,y), \quad y(x_0 ) = y_0 \) and its numerical approximation yn at the point x = xn is given by

\[ E_n = \phi (x_n ) - y_n \]
is known as the global truncation error. It arises entirely from the first two error sources listed above---that is, by applying an approximate formula to approximate data.

However, in reality we must carry out the calculations using finite precision arithmetic because of the limited accuracy of any computing equipment. This means that we can keep only a finite number of digits at each step. This leads to a round-off error Rn defined by

\[ R_n = y_n - Y_n , \]
where Yn is the value actually computed from the given numerical method. The absolute value of the total error in computing \( \phi (x_n ) \) is given by
\begin{align*} \left\vert \phi (x_n ) - Y_n \right\vert &= \left\vert \phi (x_n ) -y_n + y_n - Y_n \right\vert \le \left\vert \phi (x_n ) - y_n \right\vert + \left\vert y_n - Y_n \right\vert \\ &\le \left\vert E_n \right\vert + \left\vert R_n \right\vert . \end{align*}

Therefore, the total error is bounded by the sum of the absolute values of the global truncation and round-off errors. The round-off error is very difficult to analyze and it is beyond of the scope of our topic. To estimate the global error, we need one more definition.

As we carry the solution over many steps we would expect the values yn and \( \phi (x_n ) \) to diverge further and further apart. The local truncation error is concerned only with the direvgence produced within the present step so that it is appropriate to reinitialize yn to the value of \( \phi (x_{n} ) \) in studying this source of error. The corresponding solution to the initial value problem

\[ u' = f(x,u) , \qquad u(x_{n} ) = y_{n} \]
on the mesh interval \( [x_{n} , x_{n+1} ] \) is called the reference solution. So the difference \( u( x_{n+1} ) - y_{n+1} = h\, T_n (h) \) is called the local truncation error. We say that the truncation error is of order p > 0 if \( T_n (h) = O \left( h^p \right) , \) that is, \( \left\vert T_n (h) \right\vert \le K\, h^p \) for some positive constant K (not depending on h but on the slope function).

 

Roundoff Errors

Floating-point numbers are omnipresent: they are used in the overwhelming majority of computations in applied mathematics and computer science.

Rounding errors are not random and are typically correlated. When rounding error leads to trouble, instead of thinking about a mysterious rounding-error accumulation, you should analyze operations that may cause this error.

IV. Regular Plotting vs Log-Plotting

Plotting the absolute value of the errors for Taylor approximations, such as the ones from the previous example, is fairly easy to do in normal plotting.

The normal plotting procedure involves taking all of the information regarding the Taylor approximations and putting them into one file, then calculate the exact solution, then find the difference between each approximation and the exact solution, and plotting the solution.

Looking at the code, the syntax is pretty much self-explanatory. The only part that may confuse students is the syntax for the plot command. The plot command need not be in the same syntax as the file has it. If you wish to, you can alter the syntax to form a curve instead of dots.

Log-plotting is not much different than plotting normally. In other programs such as Mathematica, there are a couple of fixes to make, but in Maple only one line of code is required, as seen in the file. Make sure that before you use the logplot command to open the plots package.

 

 

Now, let’s compare the exact value and the values determined by the polynomial approximations and the Runge-Kutta method.

x values Exact First Order Polynomial Second Order Polynomial Third Order Polynomial Runge-Kutta Method
0.1 1.049370088 1.050000000 1.049375000 1.049369792 1.049370096
0.2 1.097463112 1.098780488 1.097472078 1.097462488 1.097463129
0.3 1.144258727 1.146316720 1.144270767 1.144257750 1.144258752
0.4 1.189743990 1.192590526 1.189758039 1.189742647 1.189744023
0.5 1.233913263 1.237590400 1.233928208 1.233911548 1.233913304
0.6 1.276767932 1.281311430 1.276782652 1.276765849 1.276767980
0.7 1.318315972 1.323755068 1.318329371 1.318313532 1.318316028
0.8 1.358571394 1.364928769 1.358582430 1.358568614 1.358571457
0.9 1.397553600 1.404845524 1.397561307 1.397550500 1.397553670
1.0 1.435286691 1.443523310 1.435290196 1.435283295 1.435286767
Accuracy N/A 99.4261% 99.99975% 99.99976% 99.99999%

You will notice that compared to the Euler methods, these methods of approximation are much more accurate because they contains much more iterations of calculations than the Euler methods, which increases the accuracy of the resulting y-value for each of the above methods.

 

  1. Gezerlis, A, Williams, M., Six-textbook mistakes in computational physics, American Journal of Physics, 2021, Vol. 89, No. 1, pp. 51--60.
  2. Higham, N.J., Accuracy and Stability of Numerical Algorithms, 2002, 2nd edition, SIAM (Society for Industrial and Applied Mathematics), Pennsylvania.

 

Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)