American Mathematical Society

The Propagation of Error in Numerical Integrations
Author(s): Mark Lotkin

Source: Proceedings of the American Mathematical Society, Vol. 5, No. 6 (Dec., 1954), pp. 869-
887
Published by: American Mathematical Society
Stable URL: http://www.jstor.org/stable/2032552 .
Accessed: 07/10/2014 04:34
Your use of the JSTOR archive indicates your acceptance of the Terms & Conditions of Use, available at .
http://www.jstor.org/page/info/about/policies/terms.jsp
.
JSTOR is a not-for-profit service that helps scholars, researchers, and students discover, use, and build upon a wide range of
content in a trusted digital archive. We use information technology and tools to increase productivity and facilitate new forms
of scholarship. For more information about JSTOR, please contact support@jstor.org.
.
American Mathematical Society is collaborating with JSTOR to digitize, preserve and extend access to
Proceedings of the American Mathematical Society.
http://www.jstor.org
This content downloaded from 193.205.210.43 on Tue, 7 Oct 2014 04:34:05 AM
All use subject to JSTOR Terms and Conditions
THE PROPAGATION OF ERROR IN NUMERICAL
INTEGRATIONS
MARK LOTKIN
1. Introduction. The numerical integration of differential equations
is generally performed by replacing the differential equations by
approximate difference equations whose solutions are expected to ap-
proach those of the associated differential equations as the step size
approaches zero. The replacement of differential by difference equa-
tions may clearly be carried out in a variety of ways; the actual
choice will depend on particular circumstances, accuracy require-
ments, computational facilities, etc.
It is now a well known fact that whenever the order of the differ-
ence equations exceeds that of the original differential equations there
are introduced certain numerical solutions that are extraneous to the
original differential equations. The behavior of these extraneous
solutions in general determines the usefulness of the integration
method. For such a method to be effective it must be "stable" in the
sense that the extraneous solutions always remain of negligible size
as compared with the actual solutions.
Thus it is of interest to distinguish first between stable and un-
stable methods. In addition, it is also of interest to determine, in either
case, the growth of error in the large, since the knowledge of this
quantity permits an estimation of the accuracy obtained.
This paper, then, deals with a number of standard methods of
integration, and investigates their stability and propagation of error.
Round-off is considered to a certain extent, but not completely; see
footnote 1. While some of the results have been obtained previously,
mainly by L. H. Thomas [1] and H. Rutishauser [2], others do not
seem to be as well known.
The propagation of error was already treated previously by
Rademacher [3
]
and others. While the method of adjoint differential
equations employed there seems to be capable of general application,
it was used, in [3] especially, for Heun's method only.
Finally there are carried out a few illustrative examples; they show
that the theoretical expressions obtained frequently lead to good
estimates.
2. Solutions of linear difference equations. Let us assume that the
Presented to the Society, September 1, 1953; received by the editors November
11, 1953 and, in revised form, May 3, 1954.
869
870 MARK LOTKIN [December
integration of the nth order differential equation has proceeded from
a starting point xo to a point Xk, and that the original differential
equation has been replaced by a difference equation of order s:
(2.1) aLkBVk+8 + ak,S-uVk+8-1 + * * * + aXk,OVk + (Xk
= 0,
where the coefficients akj, ak are known, with ak8PO for each k, and
initial values vo, v1, * * *,
v8,1
have been supplied. For our purposes
it is now convenient to use matrix notation. Introducing, then, the
column matrices
Vk+s-1 rk Vs-1
Vk+8-2 O V8-2
Uk -, bk = Uo=
L Vk _I, Vo _
and the square matrices of order s:
r-aks 0 ...0
ak,s1 ak,8-2 akO]
I 0
...
0
Jk A[ , Ak = 0 1 *.* * 0
_ 0
.
._ .
o 1 *... 1 0 O
we may express (2.1) as a system of difference equations of the first
order:
(2.2) JkEUk
=
Akuk + bk,
with E denoting the displacement operator
Evi
=
vi+.
Since Jk is nonsingular, (2.2) may be written as
-1 -1
(2.3)
Uk+1
=
Jk AkUk + Jk bk
= CkUk
+ dk,
with
-1 -
Ck = Jk
Ak,
dk
=Jk bk.
The general solution of (2.3) is
(2.4) Uk
=
Pk-1
0
+ Pt dt
twe
where
'9541 PROPAGATION OF ERROR IN NUMERICAL INTEGRATIONS 871
Pk
=-
CkPk-1, Po
=
Co,
or
k
Pk-1 =H7 Ck-t.
t=l
In particular, if Ck is actually independent of k, then Pk-l= Ck and
/ k-1
(2.5)
Uk
= Ck to + E
C-t-ldt)
t=O
The use of Sylvester's theorem [4] now permits us to express the
solution Uk in slightly different form, as follows: let G(X) =XI- C,
Ga(X) denote the adjoint of G(X), 3(X) the determinant of G(X), and
&'(X)
=
dl/dX.
If the characteristic roots Xrn, m- 1, 2, , s, of C are
distinct, then
(2.6) k
ZXHM
m=1
with
Hm -H(Xm)
=
Ga(Xtn)/a3(Xm)
Consequently, from (2.5),
"IC =
?mHm o + E Amt)'dt
?n=1 t=O
or, if u=O-,
(2.7) Uk k Xk dt.
rn=1 t=O
If the distinct roots
Xi,
i= 1, 2, , q,
of C have
multiplicities
pi,
then (2.6) must be replaced by
(2.8)
Ck
Z
1
F d"-1 (XG,(X) -]
-
di(X)
- (X -
X))ii(X).
In particular, if the only multiple root of C has the value
zero,
then
clearly (2.8) reverts to the form
(2.6)
with the summation to be ex-
tended over all the
nonvanishing
roots. This observation will be
put
to use in the subsequent discussion.
3. The variational difference equations. The
foregoing
treatment
of the difference equation will now be applied to the numerical inte-
gration of the nth order differential equation
(3.1) y(n)
=
fn(x,
y y' * * ,
y(n-1))
subject to suitable boundary conditions. The exact solution y(x) of
(3.1), assumed to exist uniquely, may be expressed in the form [2]
y(Xk+l)
= , aoiy(xk) + h 2 aljy'(xkj) + *
j=O j=-l
(3 . 2) + hN G
aNjy(N)(xk-J)
+ Tk,
y (Xk+1) = ? ai( y (XkJ) + h E ai+l,iy (Xk-;) +
j=0 -1
+ h
f
aNy
( (Xk- ) +
Tk;i
for i= 1,
2,2,
n n-1. Here h =xk+l-xk is the step used in the inte-
gration, the a(> are constants that will in general depend on r, where
r itself indicates a certain range of points
Xk_j
to the left of Xk; to each
numerical method of integration there is associated a fixed value of r.
The functions Tk, Tki are truncation errors of orders hN+1, hN-i+l, re-
spectively, with N>n denoting a positive integer. In case N exceeds
the order n of the equation (3.1) the derivatives of orders n+l,
n+2, * * -, N occurring in (3.2) may be obtained by N-n successive
differentiations of (3.1):
(3.3) y(i)(x)
=
fi(X, y(X), y'(X), * y(i-)(X)),
i=n+1, , N. The coefficients a( are not arbitrary but normally
depend on certain conditions (4.5) derived below.
Insolving the integration problem numerically (3.1) is actually
replaced by
*y
(n) =
*fn(Xf *yf *yf *y. n-1y))
where the asterisks indicate rounded values, and the method of inte-
gration actually employed may be of the form
*) *) ()
N-i
'
(5) *
(N)
(3.2a) *yk+1
E
Laj
D
*Yk-j
D .. .
i
E
*aN (0 Yk-j
j=0 j=-1
for i=O, 1, 2, , n-1, and
(3 . 3a)
*Yk+= *fi(xk+l, *yk+l... *yk+l) n
?
i
<
N.
'9541 PROPAGATION OF ERROR IN NUMERICAL INTEGRATIONS 873
In (3.2a) the circled symbols indicate pseudo-addition and pseudo-
multiplication, i.e. certain digital operations by which the cor-
responding arithmetical operations must be replaced whenever nu-
merical calculations (which of necessity involve rounding) are carried
out.'
In practice the numerical solution of the sets (3.2a) and (3.3a) is
obtained iteratively in the following manner: Extrapolation, or some
other means, permits the determination of a first set of values for
*Yk+Q,
n_i?N. Then (3.2a), with i=n-1, n-2, 0, , lead to a
first set of values for *y QF1, 0<i<? n-1. Next an improved set
*yk+,
n
<
i!<N, is computed by means of (3.3a), etc., this cycle being
repeated until duplication occurs.
Our main interest is now the determination of the errors
Ui)
*
-i
(i)
77P =yp- y ( x),
and of the associated property of "numerical stability" in the sense
that all the
77(f)v
remain small throughout the entire region of integra-
tion. By (3.2a) and (3.2),
) E [*afi) (D -
aii
y
(Xk-j)
+ * Tk,.
i=o
However, *a O *y =a*y+(*aO3*y-*a*y) +(*a-a)*y, and *aO *y
-*a*y=p, *a-a==a, with
p
P|< ug, 1 1<,u u=2-'-3Y denoting the
basic rounding error of a computation carried out to y places in a
number system of base 3.
Consequently
(.)
r
(i (X)
r
(i) (Z+ 1)
flk+l - E ai2 flk-j + ai
E
,
a, j?1fk- +
N-i
r
(i) (N)
3 . 4) +
h
aNjlk-j
-
Tki + Tki,
(3.4)+IN$f
)()
Tki
=
>(Piik + ?i j(Yk-i)v
Further, by (3.3) and (3.3a),
1
Note that there is no provision in formula (3.2a) for rounding the product in-
volving hN-i. Thus the formula and later special cases of it in ?4 are based on the
assumption that the independent variable x and the step size h may be chosen exactly,
and that multiplication by powers of h does not necessitate rounding. These assump-
tions could be dropped at the expense of including additional rounding items. As
written, however, the conclusions of the paper may not be precisely applicable to
numerical integrations in which the term in question is rounded.
7k+l= *fi(Xk+l *yk+l,
. . .
fi(xk+l, y(Xk+l), )
for i
=
n, n + 1,** ,N.
But
*fi(x, *y,
) = f(x, *y, )
+ ci,
with qi denoting certain quantities that depend on the procedure
employed in calculating
fi
from its arguments x, *y, . Thus,
correct to terms of the first order in 7,
(3.5) n7+l= (afl/ay )flk+l + 'i,
n < i <
N,
the partial derivatives to be evaluated at Xk+1, *yk+1,
. .
.,
Y(+1
The system (3.4) and (3.5) of difference equations is transformed
into the form (2.2) by the introduction of the column matrix Uk,
where
T (N) (N)
Uk = [ ?fk, 7fk-1, * * 7k-r; 1k, * * 7k-r; * 7k
,
* * 7k-r],
the superscript T denoting transposition. Then our system becomes
(3.6) JkEUk
=
AUk + bk,
where Jk, A are square matrices of order s=(r+l)(N+1), and bk
is a column matrix of the same order. These matrices are composed
as follows:
(i) The elements
Jij
of Jk are square matrices of order r+1,
Jii=I for
i=O, 1, * * *, N,
Jij=O for O<i<n-1,j<i, and for n<i<N,j>i,
= ] = for O
<
i _ n-1, 1
<
j
_ N, i < j,
L_ O
.
O_.--
-ofi/'ayi
O *
..O
0*
Jii = forn
?
i< N,O
<
j<
N-1,
i> j;
_ O * -O
(ii)
the elements A
ij
of A are square matrices of order r+1,
'9541
PROPAGATION OF ERROR IN NUMERICAL INTEGRATIONS 875
(i) (i) (i)
aio ail
...
air
1 0 0
Aii
0 0 for 0 ? i
<
n-1,
_0 0 1 0
_j
1 0 ...0
Aii
=0 for n
<
i
<
N,
_0 0 1 0_
(i) (i) (i)
ajo ail
ai,
0 0
Ais= h~i for
0
<
i
<
-1,
1
<
j
<
N,
i <
j,
LO *..O
_
Ai1
= 0 for n
<
i
<
N, j 7 i, and also for 0 < i
<
n-1, j < i;
(iii) the column matrix bk has the transpose
T
bk =
-T7, +
7r, 0,
*
0, ***
-Tk,n-1
+ 0T,n1 I0** ; f ny 01 * * 0; * ON * o} f I}
01*
.?
Clearly the determinant D(J) of J is of the form
D(J)
=
1 + Clh + C2
h2
+ **
i.e. D(J) XO for sufficiently small h.
Let us now assume that the partial derivatives
(9fi/ay(i))
have the
property that there exist points ao,
O(J)
in a suitably defined space
Ix-xoj ?a,
Iy)(-
y(J)I
<
0(i),
j=O, 1, * * , N, such that
ofil/oy(i
is
approximately equal to a constant
fij
(ao, i%,
f,
* * 1)). In
such a case let Jo denote the matrix J whose elements are evaluated
at ao,
0(j).
If the characteristic roots X of Co
=
JO 'A are distinct, then
by (2.7)
k-i
k
-
1
(3.7) Uk
=
Z Hm XM
n
o J0bt.
m t=O
Now A has at least N - n +
1 rows of zeros. There are thus
(N-n + 1) (r + 1) artificial characteristic roots of CO
of value zero.
There must be, further, the n roots associated with the n independent
solutions of the variational equations; these roots are of the form
XM= 1 +7mIh + , m = 1, 2, , n.
Thus there remain
(r + 1)(N + 1)- (r + 1)(N-n + 1)-n
=
nr
ad'ditional roots. These are "extraneous," introduced by the method
of integration. "Extraneous" solutions of the integration problem are,
consequently, solutions belonging to extraneous roots of
CO.
Now in the solution vector Uk the only components of interest are
those of 9k, where
T w ~~~~(N)
Qk =
[qk, 'q k,
. . .
X 7k ]
These may be calculated obviously from (3.7) by simply replacing
bk in (3.7) by Ck where
Ck =
[-
T + r,
-
Tkl + 7rkl, Tk,n-1 + 7rk,n-l; O)n, * *XN],
accompanied by a similar contraction of the matrices H(X) and JO1.
The determination of the characteristic values )X of CO
= JO 'A and
the construction of the error vector may be simplified somewhat,
as follows:
Let us define
(3.8) r(A) = Al - JaA,
where J, A are square matrices of orders s, and Ja denotes the adjoint
of J. Similarly, let Ga denote the adjoint of G(X) =XI-J-1A. Let,
further,
A(A)
=
det r(A),
and, as before,
8
(X)
= det G (X).
Then we have the following
LEMMA.
Ga/l'
=
ra/A'.
PROOF. Clearly
(3.9) r1(A)
=
D [AD-'
-
J-'A] = DG(A/D),
where again D = det J, G (X)
=X I- J-'A. Consequently,
A(A) = D85(A/D).
1954] PROPAGATION OF ERROR IN NUMERICAL INTEGRATIONS 877
To each root A of A (A) = 0 there is then associated a root
X = A/D
of (X) =O0. Further,
(3.10)
However, by (3.9),
(3.11) Pa(A) = D8lGa(X)
which, together with (3.9), proves the lemma.
We have thus obtained the following general
PROPAGATION THEOREM.
(3.12)
J,=
>
[Pa(A)Ja/DA'(A,)
]Z
E > Ct.
m A
In this theorem the index m is to be extended over all distinct non-
zero characteristic roots Am of A (A) =0, Xm= Am/D(J), and the ele-
ments of the vector Ck are due to truncation error and rounding. The
theorem shows again that for a method to be stable for sufficiently
small h it is sufficient that all characteristic roots Xm be of absolute
value not exceeding unity.
The actual computation of Qk would then proceed in obvious
fashion from the construction of JaA and A(A) to the calculation of
D(J), Am, and Pa(tAm) Ja, and could be carried out concurrently with
the integration.
4. The propagation of error in the case n =-N= 1. The deductions
of the previous sections will be applied now to a number of well
known methods of numerical integration. We shall start by consider-
ing the general first order differential equation
Y'
=
f(X, y).
In this case the associated homogeneous variational equation is
'7
=
AUn;
it has the fundamental solution
,q(x)
=
exp
(fyd)
Among the solutions X of the characteristic equation 8(X) =0 in-
herent in any useful method of integration there must be
one,X=Xl,
that approaches this fundamental solution 71(x) as the step size h
goes to zero. It will be seen that this root '1 is of the form
X1 = 1 +
hJvf + * exp hf-exp p, p=
hf,,;
it gives rise in (3.12) to the principal term of the form
Qk(l)
=
Mi1x Ml exp(
ffdt).
xt+1
In order to prevent the error in the large from increasing rapidly
it is thus sufficient to carry out the integration in the direction A\x
in which
f,Ax
is nonpositive.
Let us consider first the important case n =N= 1. Here we have
A Ao
h
[0 I
0
...
0]
=~Io 1ioii /~T
--JO *-
Daoo ao-a- aor
-Aoo Aol-
I 0 ..O
L K A J011 __
ho K
O 0 1 0
haio han
..
hair-
-0 ..* O0-
O *I*0 1 0 0
Ao,= All=
_0 * *
. .
_._ 1 0_
It follows that
D
=
1-Pa-,
a,ai,,-,,
1
.
0 * A0
[Jlo K]
(4.1)
0 0 ** D
r
(A)
=
Al1 JaA = [r?
]Pl]
X
~A
aoo -ao,
*
-ao,,r-i
-aOr-
-D At * O
r'oo
=
'9541
-haio -hall -hair]
(4.1) = - L a1
D0 0 0
The = LAi-alo -pai.*.*. -pai,r-i -palTj
_ 0 . -D
. .
Techaracteristic equation A\(A) =0O may be expressed in the form
(4.2) A(A) = Ar+lA0(A) =
0
A -() - AiT+l
-
eo
alr r-1 -
where
(4.3) ep= DP(ao + pa1), p = 0, 1, , r.
The equation for the extraneous roots Atm, if any, is now quite
easily obtained from (4.2).
It is convenient to write the characteristic
expression
in the form
42c
=(A)
=
Ar+ - AcO(A) - c
where
(4.4)
Aep(A)
= A DaPO +
ArolDai
+
p+
Drair,
=
= 0, 1.
Since A1=DX1-(1+p)D is a solution of A,(A) =0, there are ob-
tained for the coefficients am1 the relationships
Z aor
-
1,
(4.5)
a4 +
Z
i =- a
(j)aol = 1.
One may show, after some lengthy
calculations,
that the con-
tracted adjoint of fe(A) may be expressed as
ra(A) = AF Aco(A) hAci(A)
L f4Aco(A)
PAci(A)J
or
ra(A)
= Ar[ ]
[AcO Aci] [ j
Therefore,
ra(A)Ja
=
Ar
[Ar+l,
h(alAco(A)
+
Ac1(A))
1,
and, finally,
(4.6)
F1
]
1
[Ain, h(al\eo + Acl)]
Z
c.
Lf7ik
D m AmA' (Am)
The error
?q
may
thus be obtained
quite simply
from -lk by multi-
plication with f,.
The contribution of the root
A1,=
(1 +p)D to the error tk may now
be written down at once; it is
f7k(Al)
=Al)
[1
+ p(r
-
(r
-
l)al),
h
a,
+
E
at
l
* exp( ffdx)
c
It is of interest to apply above deductions to some specific cases.
I. Case r =0. As was pointed out above no extraneous solutions
arise in this case, so that the methods are stable in the direction in
which p <0. Techniques falling into this class are due to Euler,
Heun, Runge-Kutta, Milne, and others.
Since now
vc(A)
=
A
-
(aoo + paio),
there is obtained from (4.7) and (4.5)
k = [1 + pal, h] E exp ( f
vdx) Ct
[-
Tt
+,t]
Ct
=
or
(4.8) tlk [(-Tt + rg)/h +
alrjfv
+ cp] exp(
fydx)
dt.
xo
t?1
1, 1. Euler's method.
Yk+i = yk +
hMyk,
(aoo a1
alo)
= (1, 0, 1), T = (h2/2)y".
flk ~f
[-(h/2)y7
+ Tg/h +
i] exp( fvdx)
di.
I, 2. Heun's (modified Euler) method.
yk+l
=
Yk + (h/2)(*yk + *y'+),
(aoo a, aio)
=
(1, 1/2, 1/2),
T = - (h3/12)y"'.
b
k [-1 +
p/2, h] exp
( ffdx) ct.
II. Case r
=
1. There is one extraneous root A2; it satisfies the
equation
(4.9) Ac(A)
=
A2-Aco-pAcl
= O
where
Aco((A)
= Aaoo + Daol, Ae1(A)
= Aa1o + Dall.
Since
A, +A2=aoo+Pajo,
and ai1 satisfy (4.5), it follows that
(4.10) A2
=-
aoi
+ p(aoi
-
all).
Further
Ac(Al)
=
(2
-
aoo) + p[2(1
-
a,)
-
alo], AC(A2)
=
-A(Al).
II, 1. Simple central difference method.
Yk+1
=
yk-1 + 2lh*yk,
(aoo aoi a1 aio a1l)
= (0, 1, 0, 2, O),
T
(h13/3)y"'.
Thus
D(J) = 1,
Am=
+ I
+ P
m = 1, 2,
Ac(A1)
=
2,
and, consequently,
k= { [1 + p, 2h] E exp (f adx)
(4.11)
2
+
[-
1+ P,
2h] (-1)k-t exp( fhdx)} c.
The extraneous root A2 may thus give rise to an oscillating term
of increasing magnitude whenever the integration is applied in a
direction in which
hf,
<0. However, it is entirely possible that the
actual increase of this tierm is choked off by the rounding procedure
itself. See footnote 1. Used for integrations in the opposite direction,
over short ranges, the method may give useful results.
II, 2. Simpson's method. Here
Yk+1
=
yk-1
+ (h/3)
(*yk+
+
4*y'
+
*y'
),
whence
(aoo
aoi a1
alo a1l)
=
(0, 1, 1/3, 4/3, 1/3),
T
=
-
(15/90)yV.
It follows that
D(J)
=
1 -
p13,
Am = ? 1 + 2p/3, A (A1)
= 2.
Therefore,
k I{[1 + p, 2h] exp f Adx)
(4.12)
+ [-1 + p/3, 2h/3] (_1)k-t
*exp (- f fdx/3)} Ct.
The second root A2 may thus again make this integration method
unsuitable.
III. Case r=2. The three roots
Am,
m=1, 2, 3, satisfy the char-
acteristic equation
Ac(A)
As - Ao(A) -
PAc1(A)
= 0,
with
Aci (A)
=
A2aio
+
ADail
+
D2ai2,
i = 0, 1.
'9541
The 2 extraneous roots A2, A3 may thus be obtained from
A2 - A[(ao0
- 1) + P(a, + alo
- 1)]
-D[- aO2 +
P(ao 2- a2)] 0-.
III, 1. Adams method.
1k+l= yk + hk[*y' + (1/2)V*y' + (5/12)V
"*
y] +
with
V(i+l)qk
=
V(i)qk -
V()qk_-*
Thus, alternately,
Y*+i
yk + (h/12) [23*y'
-
16*yi
?
5*yk-2] +
Then,
(aoo, aol, ao2, a,, alo, a1l, a12)
=
(1/12)(12, 0, 0, 0, 23, -16, 5).
The 2 extraneous roots A2, A3 satisfy
A2
-(11p/12)A + (5p/12)
=
0.
Therefore,
Am = +
(- 5p/12)1"2 + llp/24,
m =
2, 3.
For sufficiently small h this method, then, is stable, and the propa-
gation of error depends essentially on A1.
Since
A
(Al)=
1 + 3p/2,
it is found that
(4.14) Xk
=
[1 + p/2, h]Zexp(ffdx )ct.
III, 2. Gregory's method. Here
yk
=
*yk + ?t
[*yk?l
-
(1/2)V'*yk?l
-
(1/12)V"*y?k+
-
(1/24)V
.
*y'i] +
-
yk + (h/24) [9*yf+ +
19*y'
-
5*yk-' +
Yk-2I
+?
.
Thus
(aoo
ao, ao2 a1 aio a11 a12)
=
(1/24)(24, 0, 0, 9, 19,
-5, 1),
T - - (19/720)h5yv.
By (4.13), then, D = 1-3p/8,
A2 - (p/6)A + p/24 = 0,
so that
Am = ?
(-p/24)112 + p/12,
m = 2, 3.
The method thus has the same stability properties as Adams'
method.
Since again
AI(Al)
= 1 + 3p/2,
one obtains from (4.7)
(4.15) 77k
=
[I
+
p/8, h] Eexp(
fAdx)
ct.
5. Numerical example. To test the propagation theorem let us
integrate the differential equation
(5.1) yI =y-2x/y
by means of Simpson's method II, 2. Taking h = 0.5 and starting at
x=O with values computed from the exact solution
y(x)
=
(2x
+
1)1/2,
there are obtained the "solutions" shown in columns (2) and (3) of
Table I. At each step a sufficient number of iterations is carried out
in order to achieve agreement to five decimals (column (2)), or four
decimals (column (3)). Due to the instability of the method the
five-decimal "solution" diverges more and more from the four-
decimal "solution."
The fourth column (4) contains the exact solution y(x), and the
fifth column (5) the error
q=*y-y(x),
the solution *y taken from
column (3).
The growth of error may be inferred from (4.12), or, somewhat
more accurately, from
?lk = k(Xl) + fk(X2),
1 k-i
-
(5.2) 27k(Xl)
-
E [(1
+
p)((hz/90)y"
+
Tt)
+
2hq]l X
2
t-=o
1
k-1
27lk(X2)
- -
E [(-1 + p/3)
((I1/90)Ay
+ rt)
(5.3)
2
(-o
h
+ (2/3)Izc]X2"1
TABLE I. INTEGRATION OF
y' =y
- 2x/y, h=0.5
(1)
(2) (3) (4) (5)
x y y
y(x) X
0 1.00000 1.0000 1.0000 0
.5 1.41421 1.4142 1.4142 0
1.0 1.73516 1.7352 1.7320 .0032
1.5 2.00529 2.0053 2.0000 .0053
2.0 2.25064 2.2507 2.2361 .0146
2.5 2.48438 2.4845 2.4495 .0350
3.0 2.73397 2.7342 2.6458 .0884
3.5 3.04775 3.0483 2.8284 .220
4.0 3.53706 3.5383 3.0000 .538
4.5 4.42054 4.4232 3.1623 1.26
5.0 6.07815 6.0834 3.3166 2.77
5.5 9.08576 9.0953 3.4641 5.63
6.0 14.31274 14.3292 3.6056 10.7
6.5 23.14506 23.1727 3.7417 19.4
7.0 37.86292 37.9089 3.8730 34.0
7.5 62.23708 62.3131 4.0000 58.3
8.0 102.49977 102.6253 4.1231 98.5
8.5 168.93852 169.1431 4.2426 165
9.0 278.52584 278.8552 4.3589 274
9.5 459.25450 459.8020 4.4721 455
10.0 757.28847 758.1877 4.5826 754
For our example,
k
= 20, h =
1/2, p5 hf = 1- [2(2x +
1)]-1,
yV
=
105(2x + 1)-91/2
= c10-5, C1
<
10,
and
X
= 1 + p.
Now p increases from 0.5 at x=0 to 0.98 at x=20, so that p-0.7
could be taken as average value. Furthermore, for the terms cor-
responding to the low values of x, which contribute most to nk(X1),
r and 4 are negligible. Thus, by (5.2),
1
19
20- t
(5.4) ?7k(Xl)
-
-
E (1/2880)yt(l. 7)
2
t-o
One obtains
7k(Xi)
-
762.
Since X2 -1 +p/3, so that
1
19
20-t
?lk(X2) t -
-
E (1/2880)yg(-O. 77) X
2 t=O
and consequently negligible, we get
11k 'qk(Xl)
-
762,
which compares very favorable indeed with the actual value of 754
for the total error.
In order to examine the oscillating term nk(X2), let us integrate
again (5.1), this time starting at x=60, and taking h= -1. The
"solution" is shown in column (2) of Table 2.
Now,
k =34, h=-1,
p =-2 + (2x + 1)-l -2,
|
T,
I
j | =
c2- 10-,
C2
2
2=
- 1
+ p/3
- 1.66.
Since for low values of t the term
(h1/90)yv,
is less then 1.10-7, in
absolute value, there is obtained the expression
33
(5.5)
|
7k(X2)
10-6 E (-
1 +
p/3)34-t.
t=O
This formula leads to
I nk(X2)
19 . 0,
which is very close indeed to the exact error given in Table 2.
The exhibited expressions, then, lead to quite useful estimates of
the error.
REFERENCES
1. L. H. Thomas, Stability of solution of partial differential equations, Symposium
on Theoretical Compressible Flow, U. S. Naval Ordnance Laboratory, 1949, pp. 83-94.
2. H. Rutishauser, Ober die Instabilitdt von Methoden zur Integration von Differ-
entialgleichungen, ZAMP vol. 3 (1952) pp. 65-74.
3. H. Rademacher, On the accumulation of errors in processes of integration, Pro-
ceedings, Symposium on Large-Scale Digital Calculating Machines, Harvard Compu-
tation Laboratory, 1948, pp. 176-185.
4. R. A. Frazer, W.
J.
Duncan, A. R. Collar, Elementary matrices, Cambridge Uni-
versity Press, 1950, p. 78.
I954] PROPAGATION OF ERROR IN NUMERICAL INTEGRATIONS 887
TABLE 2. INTEGRATION OF y'=y
-
2x/y, h= -1
(1) (2) (3) (4)
x
*Y Y(X) 105X
60 11.00000 11.00000 0
59 10.90871 10.90871 0
58 10.81665 10.81665 0
57 10.72381 10.72381 0
56 10.63014 10.63015 -1
55 10.53566 10.53565 1
54 10.44030 10.44031 -1
53 10.34409 10.34408 1
52 10.24694 10.24695 -1
51 10.14891 10.14889 2
50 10.04984 10.04988 -4
49 9.94994 9.94987 7
48 9.84875 9.84886 -11
47 9.74698 9.74679 19
46 9.64333 9.64365 -32
45 9.53994 9.53939 55
44 9.43304 9.43398 -94
43 9.32899 9.32738 161
42 9.21679 9.21954 -275
41 9.11515 9.11043 472
40 8.99193 9.00000
- 807
39 8.90203 8.88819
1384
38 8.75130 8.77496 -2366
37 8.70087 8.66025 4062
36 8.47474 8.54400 -6926
35 8.54560 8.42615 11945
34 8.10464 8.30662 -20198
33 8.53865 8.18535 35330
32 7.47987 8.06226 -58239
31 9.00031 7.93725 106306
30 6.19044 7.81025
-161981
29 11.03789 7.68115
335674
28 3.61153 7.54983
-393830
27 19.42204 7.41620
1200584
26 -12.03925 7.28011
-1931936
RCA SERVICE COMPANY, PATRICK AIR FORCE BASE

American Mathematical Society

Uploaded by

Copyright:

Available Formats

American Mathematical Society

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

American Mathematical Society

Uploaded by

Copyright:

Available Formats

The Propagation of Error in Numerical Integrations

Author(s): Mark Lotkin

You might also like