- Email: [email protected]
Chapter 5 Applications to Numerical Analysis
CHAPTER 5
5.0.
Applications to Numerical Analysis
Introduction
Despite the fact that the two theories have been developed almost independently, the connections between numerical analysis and the theory of difference equations are several. In this chapter, we shall explore some of these connections. In Sections 5.1 to 5.3, iterative methods for solving nonlinear equations are discussed and the importance of employing the theory of difference inequalities is stressed. Sections 5.4 and 5.5 deal with certain classical algorithms from the point of view of the theory of difference equations. Sections 5.6 and 5.7 are devoted to the study of monotone iterative techniques, which offer monotone sequences that converge to multiple solutions of nonlinear equations. This study also includes extension of monotone iterative methods to nonlinear equations with singular linear part as well as applications to numerical analysis. In Section 5.8 we provide problems of interest.
5.1.
Iterative Methods
By iterative methods one usually means methods by which one is able to reach a root of a linear or nonlinear equation. Let us consider the problem of solving the algebraic equation of the form F(x) = 0,
(5.1.1)
127
128
5.
where F: R
S
~
R
S ,
Applications to Numerical Analysis
which is usually transformed into the iterative form:
Xno = Xo. (5.1.2) The function f(x) is called the iterative function and is defined such that the fixed points of f are solutions of (5.1.1). There are many ways to define the iterative function f. The essential criterion is that the root x* of (5.1.1), which we are interested in, is asymptotically stable for (5.1.2), although it may not be sufficient to ensure rapidity of convergence. This, however, will imply that if Xo lies in the region of asymptotic stability of x*, the sequence will converge. If f is linear, we know that asymptotic stability can be recognized by looking at the eigenvalues of A (the iterative matrix) and indeed this is what is done in the study of linear iterative methods such as Jacobi, Gauss-Seidel, SOR, etc. Moreover, in the linear case, asymptotic stability implies global asymptotic stability. For nonlinear equations the situation becomes more difficult. In this case there are essentially three kinds of results that one can discuss: Xn+l
= f(x n ) ,
A. LOCAL RESULTS. These results ensure asymptotic stability of x*, but they have nothing to say on the region of asymptotic stability. Usually theorems in this category start with the unpleasant expression "If Xo is sufficiently near to x* ...." B. SEMI LOCAL RESULTS. The results in this category verify that an auxiliary positive function (usually the norm of the difference of two consecutive terms of the sequence) is decreasing along the sequence itself. The sequence is supposed to (or one proves that it) lie in a closed set D c R Then one can infer that the auxiliary function has a minimum in D and this minimum is located at x*. We shall see that this usually requires that x* is exponentially stable. Connected to these types of results are those requiring the stronger condition of contractivity off in a closed set D c R Contractivity implies exponential stability of the fixed point x*. C. GLOBAL RESULTS. These results say that the sequence x, given by (5.1.2) is almost (except, on a set of measure zero) globally convergent. S
•
S
•
Of course, results in the class (C) are very few in the nonlinear case. One such result is due to Barna and Smale, which says that for a polynomial with only real root, Newton's method is almost globally convergent, which we shall not present here. We shall, however, discuss in the next two sections some of the most important results in the classes (A) and (B).
5.2. Local Results
"
Let us begin with the following main result. Here f: R x* = f(x*).
S
~
R and x* satisfies S
5.2.
129
Local Results
Theorem 5.2.1. Suppose that f is differentiable with Lipshitz derivative, and the spectral radius off'(x*) is less than one. Then x* is asymptotically stable. Proof. By the transformation eTj difference equation en + 1
= f(x n) - f(x*) = f'(x*)en + =
where
f'(x*)en
=
xTj - x*, from (5.1.2) we obtain the
L
{f'(x*
+ s(xn -
x") - f'(x*)} ds en
(5.2.1)
+ g(en)
l~~ II~~:II)II = 0,
(see problem 5.2). The equation (5.2.1) is in the
form required by Corollary 4.7.2, which assures the exponential asymptotic stability of x*. Note that the exponential asymptotic stability does not imply in general that /lxl - x*11 < Ilxo - x"], unless (see section 4.2) aT] < 1, where the constant a depends on the norm used. The result established in Chapter 4 can be used to prove the next theorem on the convergence of nonstationary iterative methods defined by (5.2.2) where q : N;o x D c R ~ R , XO E D. It is obvious that nonstationary iterative methods are, in our terminology, nonautonomous difference equations. • S
Theorem 5.2.2.
S
Consider the difference equation
yno = X o ,
(5.2.3)
where f: D ~ R S , and f is locally Lipshitzian. Suppose that y* E D is an asymptotic stable solution for (5.2.3) and that the solutions y(n, no, x o) of (5.2.3) and x(n, no, xo) of (5.2.2) are in D for n E N;o' Assume further that
Ilq(n, x) - f(x)11 ~ Lnllx where L; ~ 0 for n ~ suitably chosen.
Proof.
y*ll,
(5.2.4)
Then the solution x(n, no, xo) ~ y*, when
00.
Xo
is
Rewrite (5.2.2) as Xn+I
= f(x n ) + R(n, x,)
(5.2.5)
and consider (5.2.5) as the perturbed equation 4felative to (5.2.3) with (5.2.6)
130
5.
Applications to Numerical Analysis
as the perturbation, which tends to zero as x, 4.11.1. •
~
y*. Then apply Corollary
Other results of local type can be obtained by using the comparison principle given in section 1.6. For example: Theorem 5.2.3.
Let f: D c R S ~ R S continuous, g: R+ ~ R+. Suppose that
(1) there exists x*
E
Int D such that for all x
E
B(x*, ll) < D, II
> 0,
(5.2.7) Ilf(x) - x*11 :5 [x - x*11 + g(llx - x*Ii); u + g(u) is nondecreasing with respect to u and g(O) = 0;
(2) G(u) = (3) the trivial solution of
Un + 1 =
G(u n ) ,
Uo = Ilxo - x*ll,
is asymptotically stable, with U o in the domain of asymptotic stability, and u, :5 Uo for all i. Then x* is an asymptotically stable fixed point for the equation
(5.2.8) Proof. Let Xo E B(x*, ll). Then by Theorem 1.6.1, we have [x, - x*11 :5 u; and, because Un is decreasing to zero, all x, E B(x*, 8) and the sequence will tend to zero. •
The next two theorems are applications of LaSalle's Theorem 4.8.6. Here we study the convergence of secant and Newton methods on the real line. These two methods are widely used in numerical analysis to find roots of equations and their generalizations are almost countless. On the real line, these two methods are defined by (5.2.9) and (5.2.10) respectively. Iff(zk) ¥- f(Zk+I), the difference equation (5.2.9) is well defined. Suppose that f is differentiable and that a is the simple root of f( z) = O. Using the transformation (5.2.11) and the mean value theorem f(a
+ e)
= f'(a)
+ g(a,
e)e 2 ,
(5.2.12)
S.2.
Local Results
131
we write the new difference equation corresponding to (5.2.9) having the fixed point at the origin in the form (5.2.13) where _ g(a, en+l)en+ l - g(a, en)e n M( a, en, en+ l ) ( . / a + en+!) - /(a + en)
(5.2.14)
If a is a simple root and g( a, e) is continuous and bounded, then Mt;«, ek, ek-l) is continuous and bounded if ek, ek-l are small enough. The second order difference equation (5.2.13) can be transformed into a first order system by introducing the variables x, = en and Yn = en+ l • Then the system is
= Yn, Yn+l = M(a, xn,Yn)xnYn' Lyapunov function V(x, y) = Ix[q + Iylq, X n+ l
By using the
q
2':
LlV(xn, Yn) = -(1 -IM(a, x n, Yn)Ynlq)lxnl
1, one obtains
q,
which is negative if W(x,y);;;: (l-IM(a,x,y)yj) > O. Let D(1) = {(x,y)I(!xl q + Iylq)l/q < 1)}. Since y = 0 implies W(x, y) = 1 > 0, there exists some 1) > 0 such that W(x, y) 2': 0 in D( 1). If the initial values Zo and ZI are such that (x o, Yo) and (Xl> Yl) E D( 1), then (x n , Yn) will remain in D (see Corollary 4.8.5) and will converge to the maximum invariant set contained in the set E = {(x, y) E 15(1) Ix = O}. The only invariant set with Xl = 0 is the origin, so we obtain (xn, Yn) ~ (0,0) for n ~ 00. Corollary 4.8.5 can be used to study the convergence of (5.2.10). Because there is no more difficulty involved we shall consider the case where (5.2.10) is defined in R S • The root a, which we are interested in, is simple. The change of variable similar to (5.2.11) allows us to write the method in the form where M l ( en) -_ [d/(a It can be
+ en)]-l and M ( en) -_ [d/(a + en)] . 2
ax shown that for I enII <
IIMI(en)[Mz(en)en
ax
1), -
/(a
+ en)] II =:; k(1)lle~ll.
132
S.
Then, taking Vee)
= IIell, one
Applications to Numerical Analysis
can obtain
AV(e):s -(1- k( 17)llell)llell
;:: - W(llell),
showing that A V :s 0 for k( 17) II en II :s 1. Using Corollary 4.8.5, we get a region of convergence D( 170) = {x Illx - a II < 170} with 170 = mine 17, 1/ k( 17». The parameter 17 can be chosen to maximize 170 and consequently the region of convergence.
5.3. Semilocal Results Results of semilocal type are obtained either by requiring contractivity or some of its generalizations on a compact D c R or by requiring the weaker condition that some auxiliary function decreases over a sequence. We shall consider some cases of the latter type. We begin discussing the simplest of such results in some detail to bring out clearly the arguments used. S
Theorem 5.3.1. Suppose that all x, obtained by (5.1.2) lie in a compact Do c R S and that
Ilxn+2 - xn+111:s allxn+l - xnll, with a < 1. Then the sequence converges to x* Proof.
E
(5.3.1)
Do.
Setting (5.3.2)
we have, because of (5.3.1), We can now apply Theorem 1.6.1 to obtain solution of Un + 1
Yn:S Un'
where
Un
is the (5.3.4)
= aU n
provided that Yo :S Uo. Now the equation (5.3.4) has the solution u; = 0 exponentially asymptotically stable, which means (see Theorem 4.1.1) that I;:l uj is bounded. In 00
fact from (5.3.4) we get Ij=o uj =
1 I-a
--Uo.
It then follows that
Uj,j = 1,2, ... , if we choose Uo = Ilxl IIxn+p
-
Xn
II :S IIxn+p
-
Ilxj + 1 -
Xjll
:S
- xoll. Moreover, for p ~ 1, we obtain x n+p - 1 + x n+p - 1 - X n- p-2 + ... + X n+1 - X nII
p-l
:S j=O L Ilxn+j+l -
p-l
xn+J
:S j=O L un+j •
133
5.3. Semllocal Results
Because the right-hand side can be made less than an arbitrary positive quantity when n is chosen properly, this inequality shows that x, is a Cauchy sequence and thus must converge to x* in Do- From this result we can also find the distance between Xo and x*. In fact, from n-l
n-l
j=O
j=O
[x, - xoll ~ I IIXj+l - xjll ~ I
Uj>
it follows that
Ilx* -
xoll
~
1
00
I Uj j~O
=
--Uo· 1- a
(5.3.5)
Using this property one can avoid the a priori assumption that all x, lie in Do. Also, if we assume thatfis continuous, thenf(x*) = x*. Furthermore, if we assume more than (5.3.1), namely, x, y ED, Ilf(x) - f(y) II ~ k < 1, then the solution x* is unique in Do. • Using similar arguments, one can also prove the following result. Theorem 5.3.2. Suppose that xo, Xl lies in Do and Bt x«, L;'o Uj) c Do. Then the sequence converges to x* E B(xo,LJ:=o Uj)' The next theorems are on the convergence of Newton's method defined by (5.3.6) where F is the same function as in (5.1.1). The main result was first proved by Kantorovich and is usually called Newton-Kantorovich theorem. Theorem 5.3.3. ofR s,
Assume that, for X o, x, Y E D, where D is a convex subset
(l) F is differentiable in D; (2) IIF'(x) - F'(y)11 ~
(3)
II F'(xo)-lll
(4) f3'YTJ (5)
'Yllx - yll; 11F'(xo)-1 F(xo) II
~ f3,
=
TJ;
B( Xo, f31") c
D.
Then the sequence defined by (5.3.6) is well defined and converges to a root x* of F(x) = 0 contained in B(xo, I/f3'Y)· Proof.
For every
X E
Bt x«,
II {3'Y) we have
IIF'(x) - F'(xo)11 ~
'Yllx - xoll
< {3-I.
134
5.
Applications to Numerical Analysis
A bound of IIF'(x)-111 can then be obtained by using the Banach lemma (see Corollary A5.2), which yields <: IIF'(xo)-111 <: f3 x - I - IIF'(xo)-IIIIIF'(x) - F'(xo)ll- 1 - f3'Yllx - xoll'
IIF'( )-111
From this estimate and from hypothesis (4), it follows that for x E B(xo,I/f3'Y), F'(x) is nonsingular and the Newton's iterative function is defined. Moreover,
and (see problem 5.2) IIF(xn+1)11 = IIF(xn+l) - F(xn) - F'(xn)(xn+ 1 :s; hllxn+1 - x n11
2
-
xn)11
(5.3.8)
•
Substituting for F(x n + l ) in (5.3.7), we then find
IX
2
n+2 -
f3'Y Ilxn+1 - xnl1 Xn+1I :s; - 2 1 - f3'Y II.Xn+1 - Xo II'
Suppose for a moment that the sequence X n lies in B(xo, 1/ f3y), then we . gUn, ( "n "n-I) = f3'Yu~n ' The can apply Theorem 4.5.4 with L.j=1 Uj, L.j=1 U j
2(1 - f3y
Lu
j )
j=O
comparison equation is therefore f3Yu~
(5.3.9)
whose solution is Un = (,8'Y)-lr n- l. The origin is thus exponentially stable for (5.3.9). All x, lie in B(xo, 1/ ,81') because [x, - xoll < (f3'Y)-I2:;::11 2-
j
-
1
<
;1"
We may apply the
same arguments used
b ef ore. In f act, II xn+ p _ Xn II <: - "p-I L.j~O II xn+j+ 1 _ xn+j
II
<: -
"p-I _ (f3'Y)-1 L.j=O un+j 2
x 2:;::~ r n - j = (f3'Y)-lr n(l - r p ) , which shows that X n ~ x* and also gives the error estimate IIx* - xnll:s; (f3'Y)-lr n. The proof that x* is a root of F(x) follows from F(xn) = F'(Xn)(Xn+1 - x n) = [F'(xo) - (F'(x n) - F'(xO))](Xn+1 - x n) and then IIF(xn)ll:s; [IIF'(xo)11 + 'Yllxn - xoll]llxn+1 - xnll:s; [IIF'(xo)11 + 'Y(f3'Y)-I]llxn+1 - xnll. Taking the limit as n ~ 00 one has, by the continuity of F, F(x*) = O.
5.3.
Semilocal Results
135
Hence by Theorem 4.5.4, the conclusion follows. • The conclusion of Theorem 5.3.3 has been easily obtained because the solution of (5.3.9) can be given in the closed form. More refined estimates can be obtained if we assume in (5.3.9) that Uo < !(f~y)-l. In this case, the solution of (5.3.9) is more difficult to compute but it can be done. In fact, let us define the new sequence Zn by dZn- 1 = {3YUn, L l = O. Then . (5.. 3 9) IS . transtorme .. d'mto Zn+l = Zn + -1 (zn - Zn_l)2 , W h'IC h equation 2
1 - z;
possesses the first integral, that is, the sequence {z.} satisfies the first order . Zo - !z~ . of equation Zn+l = , where Zo = {3yuo . (See (1.5.12).) The solution 1 - Zn this equation is z; = 1 - (1 - 2Z0 ) 1/ 2 coth k;", (see problem 1.27) where k is a constant depending on the initial condition. (See problem 5.8 for its determination and for a more refined bound.) This solution exists if Zo s ! and then Uo s !({3y)-l. The case of equality has already been examined. Suppose that Ul < !({3y)-I. Then it follows that lim z; = 1 - (1 - 2Z0 ) 1/ 2 == z* and since for n~OO
• , we obtain z; - z* = 2(1 - 2z o) 1/2 e" k2 + . • Zn - z* Consequently, we get !~~ (Zn-l _ z*? = 1, which shows that the order of
large n, coth k2 n
= 1 + 2 e" k2
n+l
n
1
convergence of the sequence z.,, to the limit point, is quadratic. From
the
comparison
_l_(Zn - Zn-l) and hence {3y
[x, of z.;
xoll s l:;':-~ Uj S
Ilx
n -
principle
it
x*11 < ~(1 {3y
follows
that
Ilxn+1 -
Xn
II <
1
2Z0 ) 1/ 2 e- k 2 n + • Also for n > 0,
_l_ Zn_ 1 s·_l- z * because of the increasing behavior
{3y
(3y
•
Summing up these considerations we see that we have proved a variant of Theorem 5.3.3, which we state below. Theorem 5.3.3(a).
Assume that the hypotheses (1), (2), (3) of Theorem 5.3.3
hold. Suppose further that (4) f3yY]
Z*)
and (5) B( Xo, f31 c: Dare y satisfied. Then the sequence x, is well defined and converges quadratically to
x* contained in B( xo, (31y z*). Newton's method is invariant under affine transformation, that is, it is invariant under the transformation G = AF where A is an invertible matrix,
136
S. Applications to Numerical Analysis
whereas the hypotheses in the Newton-Kantarovich theorem are not. In fact, one has 0- 10 = F-1A-1AF = F- 1 F, while for example, the hypothesis (3) is not invariant under the same transformation, The following variant of Theorem 5.3.3 takes care of this situation. Theorem 5.3.3(b). and suppose that
Assume that the hypotheses (1) of Theorem 5.3.3 holds
II F'(XO)-I F(xo) II :s; a; IIF'(xo)-I(F'(y) - F'(x»11 :s; wllx - yll;
(2) (3) (4) (5)
aw:s;!; B(xo, w-1z*)
E
D, where z* = 1 - (1- 2wa)I/Z:s; 1.
Then the sequence x, remains in B(xo, w-1z*) and converges quadratically to a root x* of F(x).
Proof.
We have
Xn+Z - Xn+1 = -F'(Xn+I)-IF(xn+l)
= -F'(Xn+I)-I[F(xn) + F(xn+ 1 )
-
F(xn)]
-F'(Xn+I)-I{ F(x n) + F'(Xn)(Xn+1 - x n)
=
+ = -
f
[F'(x n + t(xn+1 - x n» - F'(x n)] dt(x n+1 - x n)}
F'(X n+1)-1 F'(xo){I F'(XO)-I[F'(xn + t(xn+1
-
»
xn
- F'(xn)] dt(x n+1 - Xn)}, from which it follows that Ilxn+z - xn+lll -s IIF'(xn+I)-1 F'(xo) II . ~ Ilxn+1 - x,
liZ,
Since F'(XO)-IF'(xn+ l ) == F'(XO)-I F'(Xn+l) - F'(xo)-I F'(xo) + F'(XO)-I F(xo), using hypothesis (3), we get 11F'(xO)-I(F'(xn+l) - F'(xo» II < wllxnn - xoll < z*:s; 1, if Xn+l E B(xo , w-1z*). Then, by the Banach lemma it follows that F'(x n+ l) is invertible and II F'(Xn+I)-1 F'(xo) II :s; 1-
II 1
W Xn + 1 -
II' Finally, we arrive at
Xo
IIXn+Z -
w Ilxn+1 - xnf xn+lll :s; -2 1 _ I _ II' w Xn+1 Xo
137
5.3. Semilocal Results
which is the same as before with f3y replaced by w. All the previous results then apply, and we are done. • The next result also shows the familiar quadratic decay of the errors in the Newton method, but it requires more assumptions on the inverse ofthe Jacobian.
Assume that F: Do c R S Moreover, suppose that
Theorem 5.3.4. (1) (2) (3) (4) (5)
~
R S and Do c R S is a convex set.
for x E Do, F is differentiable, for x E Do, F'(x) is invertible and IIF'(x)-111 for x, y E Do, 11F'(x) - F'(y)11 :5 yllx - yll, for Xo E Do, IIF'(xo)-IF(xo)II:5 TJ, and (\' = 4f3YTJ < 1.
:5
f3,
Then the Newton iterates remain in ei«; ro), where ro converge to a solution x* of F(x) = O. Proof.
= TJ 2:~~0 (\'ZH
and
From
IIx n+z -
xn+III =
IIF'(xn+ I)-IF(x n+I)11
:5
f3IIF(xn +1 ) II
and from the estimate of II F(xn + l ) I obtained in the previous theorem, we get (5.3.10) The comparison equation is therefore
f3y
Z
Un + 1 = TUn,
Uo
=
(5.3.11)
TJ,
whose solution is u; = TJ(\'zn_ l • Now applying comparison Theorem 1.6.1 we get (see problem 5.4)
[x, - xn+11 1:5
Un
= "1(\'zn_ l.
(5.3.12)
Since (\' < 1, the trivial solution 01 (5.3.11) is exponentially stable and hence it follows that x, is a Cauchy sequence. Furthermore, we also have [x, xoll :52:';:0 Uj = ro, which means that all the iterates remain in Bt x«, ro) and converge to x*. • Corollary 5.1.1.
The error [x" - Xn II satisfies the inequality
Ilx* - x, I :5 en IIxn -
Xn-I
liz,
for n
= 1,2, ....
(5.3.13)
As an application of Theorem 1.6.7, let us consider the two-step iterative method
138
5.
Applications to Numerical Analysis
Theorem 5.3.5. Suppose that G: D x D c R x R set Do c R the following inequality holds true S
S
~
R S and on some closed
S
,
IIG(x,y) - G(y,
z)ll:5 allx - yll + ~lly - zll, let there exist x o, XI E D such
where a + ~ < 1. Furthermore, that Xk E Do for k e: 1. Then the iteratives converge to the unique fixed point of G(x) = G(x, x). Proof. Setting Yk = Il xk + l - xkll, we obtain Yk+l:5 aYk + ~Yk-l' By Theorem 1.6.7, we then get Yk :5 Uk, where Uk is the solution of Uk+l = auc + ~Uk-I' which is Uk = coz~ + CIZ~, where ZI and Z2 are the roots of Z2 - az - ~ = O. Because of the assumption on a and ~, the two roots ZI and Z2 are less than 1 in modulus, which means that the zero solution of the comparison equation is exponentially stable. Using similar arguments as in the previous theorems, we can conclude that the sequence Xk is a Cauchy sequence. Since
IIG(x*) - Xk+111 :5IIG(x*, x") - G(x*, Xk-I)II + IIG(x*, Xk-I) - G(Xk' Xk-I)II :5 ~llx* - xk-t11 + allx* - xkll, the convergence of Xk shows that Xk+l ~ G(x*). The uniqueness of x* now • follows very easily and the proof is complete. The effect of perturbation on the iterative methods are of two different kinds. If the perturbations are small and tend to zero as n tends to 00, then the theorems on total stability will ensure that the fixed points continue to have asymptotic stability. In numerical analysis, however, perturbations may remain bounded for all n (roundoff errors, for example). In this case, a more convenient concept is the practical stability. If an iterative procedure is practically stable, the sequence of iterates will not tend to the solution x* but to a ball B(x*, B) surrounding x*. Inside B, the solution may oscillate but what is important is that it never leaves B. Of course it would be nice to have B as small as possible. Let us consider for example the iterative method
= Xo'
(5.3.14)
Ilf(x) - f(y) II :5 allx - yll
(5.3.15)
Xn+ 1 = f(x n ) ,
We assume that for x, y
E
Do c R
x na
S ,
with a < 1. Suppose that the errors in the computations perturb (5.3.14) by a bounded perturbation R n , that is
xn + 1 = f(x n ) + R n ,
xna = xo,
(5.3.16)
5.4.
139
Unstable Problems
with X n , xn E Do for all n. Then Theorem 4.11.2 can be applied. The conclusion is that the difference of the two solutions [x, - xn II will not exceed Rj(l- a), where R
B(Xn,~) I-a
=
sup Rn- The condition that all the balls
must be contained in Do may be very restrictive. In fact if a
is close to one, it follows that these balls are very large and the method is considered a bad one.
5.4.
Unstable Problems: Miller's, Olver's and Clenshaw's Algorithms
The situation that we shall analyze in this section and in the next is in a certain sense the opposite of that treated in the previous ones. In fact, there the limit point was important (and indeed in the theory of iterative processes the "solution" is the limit point), where the intermediate values z; are considered an unavoidable noise. Here the limit point will not be important (generally will not exist) and the important part ("the solution") will be a finite subset of the sequence itself. These problems are very typical in that part of numerical analysis which concerns the study of special functions, orthogonal polynomials, quadrature formulas, and numerical methods for ordinary differential equations. Consider the homogeneous scalar linear equation (see (2.1.5)) of order k. (5.4.1) It has k linearly independent solutions/l(n)'/in), ... '/k(n). Suppose that the solution Yn = 0 (n ~ 0) is unstable, and
. II(n) lim 1"( n ) = 0,
i = 2,3, ... , k.
(5.4.2)
n~oo.li
When this happens, the solution II(n) is said to be minimal. The problem is how to find II or a multiple of it. Even if we know the exact initial condition that generates the minimal solution, small errors in the calculations will, as usual, lead us to solve a perturbed equation LYn = en, whose solution (see (2.1.11)) will contain all the j;(n) and this will destroy the process (see for example [52], [175]), because the h(n) (i ~ 2) grow faster than II(n). The same problem appears, of course, in the nonhomogeneous case, where it is not excluded a priori that there exists a bounded solution even if all the hen) are unbounded. One is often interested in following this solution.
140
5.
Applications to Numerical Analysis
Consider, for example, the following first order equation Yn+l - (n
+ 1)Yn
= -1,
Yo = e.
(5.4.3)
The solution of the homogeneous part is Yn = en!, while the complete solution is Yn = n!( e -
2:.;=oh) ~ 2:.:=1 (n +l s
) (S) ,
which is bounded for
n ~ 00. If one tries to follow this solution, a small error in the initial data e (always existing because it has infinite digits), will be amplified by the factor n!. A simple idea is to find the way to look for the solution of the problem backwards (for n ~ -00), because in this case the origin becomes asymptotically stable for the homogeneous equation and the errors will be damped. This idea indeed works, as we will show in a simple case, which is the original proble~ to which it was applied (see [52] for references). Consider (5.4.6) where an, en ¥- 0 for all n, and the nonhomogeneous problem LYn = gn'
(5.4.7)
Suppose that two linearly independent solutions ¢n and l/Jn of the homogeneous equation are such that ¢n is minimal, that is
· ¢n 0. IIm-= "_00 t/Jn
(5.4.8)
Of course all solutions that are multiple of ¢n are minimal and vice versa. That means that only one appropriate initial condition is needed to determine the minimal solution we are interested in. The general solution of the nonhomogeneous problem is then (5.4.9) where
Yn
is a particular solution that will be supposed minimal too; that is,
· Yn 0 IIm-= I/Jn '
n_OO
(5.4.10)
and Yo = O. We are interested in the solution Yn
=
Yo,/, +_ ¢o't'n Yn'
A class of methods are based on the following result.
(5.4.11)
5.4.
141
Unstable Problems
Theorem 5.4.1.
Consider the boundary value problem Ly
= gn,
y~N)
= Yo,
Y<:')
= 0,
(5.4.12)
where N is a positive integer, and suppose that conditions (5.4.8) and (5.4.10) are verified. Then the problem (5.4.12) has a solution and moreover, for fixed n, y~N) ~ Yn as N ~ 00.
Proof. Since Yn is a particular solution of the nonhomogeneous equation, any other solution of the same equation will be y~N)
= c~N)cPn + ciN)ljJn + Yn'
(5.4.13)
The boundary conditions give for C~N) and ciN) the values YN C(N) _ I
-
IjJN
ljJo
It follows that C~N) and ci
(N) __ cPo (N)
.
cPo _ cPN'
C2
-
ljJo
CI'
IjJN
N)
tend to zero as N ~ 00 and the claim follows. Note that the condition Y<:') = 0 can be replaced by Y<:') = k where k is any fixed value (for example an approximation of YN if available). In the homogeneous case y~N) can be obtained starting from j<:') = 0, j~NlI = 1, and then using the difference equation backwards ~(N)
Yn-I
an _
b., :
c,
en
= --Yn+1 --Yn
obtaining a value y~N) different from Yo. Multiplying then the sequence by the scale factor Yolj~N), one obtains y~N). This is the Miller's algorithm. A disadvantage of this algorithm is that one does not know a priori which value of N must be used. A modified version of it, which works in the nonhomogeneous case is also obtained by considering that the boundary value problem (5.4.9) is equivalent to the system of equations Ay(N) = b, (5.4.14) where
bl
0
0
al
CI
A=
0
0 aN-I
0
0
bN - I
CN-I
and Y
(N) _
-
b
«(N)
Yt
= (gt
(N)
"'YN-t
)T
,
- cIYo, ... , gN-t) T.
142
5.
Applications to Numerical Analysis
Some clever method to solve this system can be used to avoid the growth of the errors (see [23], [24]). By solving the system in an appropriate way, it is also possible to determine the best value of N (Olver's algorithm). • Similar in motivation to the previous algorithms and indeed related to them is the following, known as Clenshaw's algorithm. Here the problem is to evaluate the sum L~=o bnun , where b; is a given sequence and u.; is supposed to satisfy a k th linear difference equation LUn = O. The algorithm uses the result stated in Theorem 2.1.8, that is, one considers the transposed equation LTYn = bn, YN+j = 0 (j = 1,2, ... , k), obtaining N
I
n=O
k-l
b.u;
=I
j=O
j
Yj
I
i=O
Pi(j - k)Uj_i'
Suppose, for example, we want to compute IN = Ln~o bnTn(z), where Tn(z) are Chebyshev polynomials satisfying the equation Tn + 1 - Zz'T; + T; = O. The transpose equation is N
Yn - 2ZYn+l YN+1
+ Yn+2
=
bn,
n
= N, N
- 1, ... , 0
= YN+2 = 0,
which can be solved recursively obtaining Yl and Yo. From the quoted result it then follows IN = Yo + Yl [T1(z) - 2zTo(z)] = Yo - ZYI' It is interesting to note that the sum can be obtained without the knowledge of the polynomials Tn(z) except for To(z) and T1(z) and also the reduction of the operations involved.
5.5.
Unstable Problems: the "SWEEP" Method
The topic discussed in this section is very similar to the previous one, but historically it has different sources. It arises typically in that part of numerical analysis that concerns the approximation of solutions of ODE and PDE. We shall consider, for simplicity, the autonomous case. The more general case is not, in principle, more complicated. Let us consider Yn+l - 25Yn + Yn-l
= gn,
y(O) = a,
Yk
= {3,
(5.5.1)
with 5 > 1, k> O. This problem can be solved in many different ways. It can be solved rewriting it as a system of k - 1 equations similar to (5.4.14). In order to control the growth of errors it is better to treat the problem as a difference equation.
5.5.
143
Unstable Problems
The characteristic polynomial of (5.5.1) has two positive roots, one less than one and the other greater than one. This means that the origin is not asymptotically stable for the homogeneous part and this, as usual, gives troubles because the errors are amplified. One way to avoid this difficulty is to transform the linear second order equation (5.5.1) into nonlinear first order equations (see problem 5.12). 1 Xn+1 = 2i> _ x ' (5.5.2) n
and (5.5.3) Furthermore, one verifies that Yn-l
= XnYn + z.:
(5.5.4)
The solution is then obtained by computing recursively the sequences {x n } , {zJ using (5.5.2) and (5.5.3) with XI = 0 and ZI = a up to n = k - 1 and then computing {Yn} using (5.5.4) from n = k to n = 1 with Yk = (3. The advantage is that (5.5.2) has a limit point asymptotically stable and the sequence x, converges monotonically to this point, while the sequences (5.5.3) and (5.5.4) remain bounded.
Theorem 5.5.1. If Xo = 0, 8> 1, the sequence {x n } converges monotonically to the first root p of the characteristic polynomial of (5.5.1) and, moreover, 0:$ x, < p. The critical points of (5.5.2) are the roots of the equation x 2 - 2i>x + 1 = 0, which is the characteristic polynomial of (5.5.1). Let p be the root less than one. Changing the variable h; = p(p - x n ) , equation (5.5.2) becomes Proof.
p2 hn + 1 = --h-hn • 1+ n
(5.5.5)
One shows at once (for example using the Lyapunov function Vn = h~) that the half line h > 0 is contained in the asymptotic stability region of the origin, and moreover 0 < h; :$ p2, from which it follows that -p :$ x, - P < O.
•
(5.5.6)
The next theorem is not really necessary in this case, but it could be useful in the more general case.
s.
144
Applications to Numerical Analysis
Theorem 5.5.2. In the hypothesis of Theorem 5.5.1 the solution x, = P is practically stable for (5.5.2). Proof.
The conditions of Theorem 4.11.2 are satisfied for x, < p.
•
The last result ensures that for equation (5.5.2) the errors will not be amplified. Let us consider the equations (5.5.3) and (5.5.4). Theorem 5.5.3. In the hypothesis of Theorem 5.5.1, the sequence {z.} remains bounded, even if the two equations are perturbed with a bounded perturbation. Proof.
From
p2 - 20p
+ 1 = 0, it follows that 20' - p = 1/ p and = P - gn( + Zn ) and f rom (5.5.6 )
. ( ) can be wntten 5.5.3 as Zn+1
IZn+tl:::; p(lgnl + IZnl)·
I - p xn - p
Using Corollary 1.6.1, it follows that IZnl:::; pna + L;:~ pn-s+1lgsl from which we see that if'[g.] is bounded z; is bounded. Clearly if we add to gn a bounded perturbation Tn (any sort of errors), the statement remains true. • Similar arguments can be applied to (5.5.4) since it follows that IYn-tl :::; plYnl + IZnl· The method just described is known as the "SWEEP" method (see [28]).
5.6. Monotone Iterative Methods In this section, we shall develop a general theory for a broad class of monotone iterations. This class of iterations includes Newton's method as well as a family of methods, which are called Newton-Gauss-Seidel processes that are obtained by using the Gauss-Seidel iteration on the linear systems of Newton's method. As before, we are interested in finding the solutions of
0= f(u), where f
E
C[R
n
;
R
n
].
(5.6.1)
Let us first split the system (5.6.1) as
0= j;(u i , [u]p" [u]q,),
(5.6.2)
where, for each i, 1:::; i:::; n, Pi + qj = n - 1 and u = (Ui, [u]p" [u]q')' Let v, W E R" be such that v :::; w. Then v, ware said to be coupled quasi lower and upper solutions of (5.6.1) if
o:::;j;(v [v]p" [w]q,), o-z.j;(w j, [w]p" [v]qJ i,
145
5.6. Monotone Iterative Methods
Coupled quasi extremal solutions of (5.6.1) can be defined easily. Vectorial inequalities mean the same inequalities between the components of the vectors. We are now in a position to prove the following result.
Theorem 5.6.1. Assume that f E C[R n , R n ] and possesses a mixed quasimonotone property, that is, for each i,};(uj, [u]p" [u]q') is non decreasing in [u]p, and nonincreasing in [u]q,. Suppose further v, ware coupled quasi lower and upper solutions of (5.6.1) and
j;(Uj, [u]p" [u]q) - };(iij, [u]p" [u]q.) 2:= -Mj(uj - iij) whenever v ~ ii -s U ~ wand M, > O. Then there exist monotone sequences {v n }, {w n } such that Vn ~ p, Wn ~ r as n ~ ex) and p, r are coupled minimal and maximal solutions of (5.6.1) such that v ~ p ~ u ~ r ~ wfor any solution u of (5.6.1). Proof.
For any
771,772 E [v, w] = [u E R" : v ~ u ~ w].
Consider the
system
u, = Mi l};( 77ti, [77I]P,' [772]q.)
+ 77ti,
i
= 1,2, ... , n.
Clearly U can be uniquely defined, given 771, 772 E [v, w]. Therefore, we can define a mapping A such that A[ 771, 772] = U. It is easy to show that A satisfies the properties: (i) 0 ~ A[v, w], 02:= A[w, v]; (ii) A is mixed monotone on [v, w]. Then the sequences {v n }, {w n } with Vo
= v, Wo = w can be defined as follows
Furthermore, it is clear that {v n } , {w n } are monotone sequences such that v ~ V n ~ W n ~ wand consequently lim V n = p, lim W n = r exist and satisfy the relations n_OO n_OO
By induction, it is also easy to show that if (Ul> U2) is a coupled quasi solution of (5.6.1) such that v ~ Ul> U2 ~ w, then V n ~ Ul> U2 ~ W n and consequently (p, r) are coupled quasi extremal solutions. Since any solution U of (5.6.1) is a coupled quasi solution, the conclusion of the theorem follows and the proof is complete. • If f does not possess the mixed quasi monotone property, we need a different set of assumptions to generate monotone sequences that converge to extremal solutions. This is the content of the following results.
146
S.
Theorem 5.6.2.
Applications to Numerical Analysis
Assume that
(i) there exist v, W E R" with v :s; W such that 0 :s; f( v), 0 =:: f( w); (ii) there is a n x n matrix M such that f(y) - f(x) =:: -M(y)(y - x), whenever v :s; x :s; y :s; w. Then the sequence {w n}, given by Wn + 1
= Wn + Bnf(wn ) ,
(5.6.3)
where B; is any nonnegative subinverse of M( wn), is well defined, and {w n} is monotone non increasing such that lim Wn = r. If there exists a nonsingular n~OO
matrix B =:: 0 such that lim inf B; =:: B, then r is a maximal solution of (5.6.1). n-e co
Set v = Vo and W = Wo' From B o > 0 and f( wo) :s; 0 it follows that Using the fact that B o is a subinverse of M (w o), we find for any [v, w],
Proof.
WI :s; W o'
u
E
u
+ Bof(u) =
WI -
:s; WI -
Hence, in particular, vo:S; Vo f( WI)
:s;
(w o - u) - Bo[f(w o) - f(u)J
[I - BoM(wo)](wo - u):s;
+ Bof( vo) :s;
f( wo) + M( woH Wo -
WI) =
WI'
WI'
Similarly, we obtain
[I - M( wo)BoJf(wo) s
o.
Proceeding similarly we see by induction that f(w n ) =:: 0,
n = 1,2, ....
(5.6.4)
Consequently, as a monotone nonincreasing sequence that is bounded below, {w n } has a limit r 2: Vo. If u is any solution of (5.6.1) such that u E [v, w J, then u = u + Bof( u) :s; WI' then by induction u :s; W n for all n. Hence u :s; r. Finally, the continuity of f and the fact lim inf B; =:: B, where B =:: 0 is nonsingular, (5.6.4) yields n
= -(lim inf B; )f( r) =:: -
Bf( r)
2:
0,
n~OO
which implies f(r) = 0 completing the proof.
•
By following the argument of Theorem 5.6.2, one can prove the next corollary.
5.7.
147
Monotone Iterative Methods
Corollary 5.6.1.
Let the assumptions of Theorem 5.6.2 hold. Suppose that M(y) is monotone non increasing in y. Then the sequence {v n } with Vo = v, given by V n+1
= Vn + Bnf(vn)
is well defined, monotone nondecreasing such that lim minimal solution of (5.6.1). n-+OO
Vn
= p, and is the
The case of most interest is when p = r, because then the sequences {v n }, {w n } constitute lower and upper bounds for the unique solution of (5.6.1). The following uniqueness result is of interest.
Theorem 5.6.3. Suppose that f(y) - f(x):5 N(x)(y - x), where N(x) is nonsingular matrix and N(X)-l ~ O. Then if (5.6.1) has either maximal or minimal solution in [v, w], then there are no other solutions in
[v, w]. Proof. Suppose r E [v, w] is the maximal solution of (5.6.1) and u is any other solution of (5.6.1). Then
0= f(r) - f(u)
:5
E [v, r]
N(u)(r - u).
Since N(U)-l ~ 0, it follows that r ~ u and hence r = u. A similar proof holds if the minimal solution exists. The proof is complete. •
5.7.
Monotone Iterative Methods (Continued)
Consider the problem of finding solutions of Ax = f(x),
(5.7.1)
where A is an n x n matrix and f E C[R n, R n ] , which arises as finite difference approximation to nonlinear differential equations. If A is nonsingular, writing F(x) = f(x) - Ax = 0, one can study existence of multiple solutions by employing the method of upper and lower solutions and the monotone iterative technique, described in Section 5.6. In this section we extend such existence results to equation (5.7.1) when A is singular. For convenience let us split the system (5.7.1) and write in the form (5.7.2)
148
5.
Applications to Numerical Analysis
where for each i, 1 ~ i ~ n, u = (Ui' [U]Pi' [u]q') with Pi + qj = n - 1 and (Au); represents the i l h component of the vector Au. As before, a function f E C[R", R"] is said to be mixed quasi-monotone if, for each i, /; is monotone non decreasing relative to [u ]p, components and monotone nonincreasing with respect to [u ]q, components. Let v, w E R" be such that v ~ w. Then v, ware said to be coupled quasi lower and upper solutions of (5.7.1) if
Coupled quasi-extremal solutions and solutions can be defined with equality holding. We are now in a position to prove the following result. Theorem 5.7.1. Assume that (i) f E C[R n , R n ] and f is mixed quasi-monotone; (ii) v, ware coupled quasi lower and upper solutions of (5.7.1); (iii) /;(u;, [u]p" [u]q') - /;(u;, [u]p" [u]q,) ~ -Mj(u j - uJ whenever v ~ u ~ u ~ wand M, > 0, for each i; (iv) A is an n x n singular matrix such that A + M = C is nonsingular where M is a diagonal matrix with M, > 0 and C- I ~ O. Then there exist monotone sequences {v n }, {w n } such that V n ~ p, W n ~ r as n ~ 00 and p, r are coupled quasi-extremal solutions of (1.1) such that if u is any solution of (1.1), then v :5 P ~ u ~ r ~ w.
Proof. For any YJ, I-t system
E
[v, w]
=
[x
E
Rn:v
~
x
~
w], consider the linear
(5.7.3)
Cu = F( YJ, I-t),
where C=A+M and for each i, Fj(YJ,I-t) =/;(YJk, [YJ]p,,[I-t]q,) + MjYJj. Clearly u can be uniquely defined given YJ, I-t E [v, w] since the matrix C is nonsingular. Furthermore, we have Av ~ F(v, w), Aw ~ F(w, v), and F is mixed monotone. Consequently we can define a mapping T such that T[ YJ, I-t] = u and show easily, using the fact that C- I ~ 0, that (a) v ~ T[v, w], w ~ T[w, v]; (b) T is mixed monotone on [v, w]. Then the sequences {v n } , {wn } with follows:
Vo
= v,
Wo
=w
can be defined as
S.7. Monotone Iterative Methods
149
Furthermore, it is evident from the properties (a) and (b) that {vn } , {w n } are monotone such that Vo ~ VI ~ ••• ~ vn ~ Wn ~ ••• ~ WI ~ Wo0 Consequently, lim e, = p, lim Wn = r exist and satisfy the relations n-JoOO
n-JoCO
By induction, it is also easy to prove that if (u l , U2) is a coupled quasisolution of (5.7.1) such that v ~ UI, U2 ~ w, then V n ~ ul , U2 ~ Wn for all nand hence (p, r) are coupled quasi-extremal solutions of (5.7.1). Since any solution u of (5.7.1) is a coupled quasi-solution, the conclusion of the theorem follows and the proof is complete. • The case of most interest is when p
=
r, because then the sequences {vn },
{w n } constitute lower and upper bounds for the unique solution of (5.7.1).
The following uniqueness result is thus of interest. Corollary 5.7.1. for each i, that
If, in addition to the hypotheses of Theorem 5.7.1, we assume (5.7.5)
where B is an n x n matrix such that (A - B) is nonsingularand (A - B)-I O. Then u = p = r is the unique solution of (5.7.1) such that v ~ u ~ w. Proof.
Since p
~
r, it is enough to prove that p
2::
2::
r. We have by (5.7.5)
[A(r - p)]j ~j;(ri' [r]pi' [p]q) - j;(Pj, [P]Pi' [r]q) ~ [B(r - p)]j.
Consequently, (A - B)(r - p)
~
0, which implies r ~ p.
•
Remark 5.7.1. If qi = 0 for each i, then Theorem 5.7.1 shows that r, pare maximal and minimal solutions of 5.7.1 and Corollary 5.7.1 gives the unique solution.
If f does not possess the mixed quasi monotone property, we need a different set of assumptions to generate monotone sequences that converge to a solution. This is the content of the next result. Theorem 5.7.2. Assume that (iv) of Theorem 5.7.1 holds. Further, suppose that V, W E R" such that v ~ W, Av
~f(v)
- B(w - v),
-B(x - y)
~f(x)
Aw 2::f(w) - f(y)
~
+ B(w -
B(x - y)
v),
(5.7.6) (5.7.7)
150
5.
Applications to Numerical Analysis
whenever v ~ y ~ x ~ w, B being an n x n matrix of nonnegative elements. Then there exist monotone sequences {v n }, {w n } that converge to a solution u of (5.7.1) such that
provided that (A - B) is a nonsingular matrix. Proof.
We define F(y, z)
= ![f(y) + fez) + B(y -
z)).
(5.7.8)
It is easy to see that F(y, z) is mixed monotone and
-B(z - z)
~
F(y, z) - F(y, z)
whenever z, z, y, Y E [v, w] and i
~
~
B(y - y)
(5.7.9)
z, Y ~ y. In particular, we have
F(y, z) - F(z, y)
= B(y -
z).
(5.7.10)
From (5.7.8) we obtain -B(w - v) ~ F(v, w) - f(v), which yields Av ~ f(v) - B(w - v) ~ F(v, w).Similarly Aw ~ F(w, v). Finally, it follows from (5.7.8) that F(x, x) = f(x). Consider now for any Tl, JL E [v, w], the linear system given by Cu = O( Tl, JL) where C = A + M, M being the diagonal matrix with M; > 0 and for each i, OJ( Tl, JL) = F; ( Tl, JL) + MjTli' Proceeding as in proof of Theorem 5.7.1, we arrive at Ap = F(p, r), Ar = F(r, p). Using (5.7.9), we see that A(r - p) = F(r, p) - F(p, r) = B(r - p) and this implies u = r = p is a solution of (5.7.1). The proof is complete. • Corollary 5.7.2. If in addition to the hypotheses of Theorem 5.7.2 we also have f(x) - fey) ~ C(y)(x - y) for x, y E [v, w], where C(y) is an n x n matrix such that [A - C(y)] is nonsingular and [A - C(y)r 1 ~ 0, then u is the unique solution of (5.7.1).
Proof.
If ii is another solution of 5.7.1, we get A(u - ii) = feu) - f(ii)
which yields u (5.7.1). •
~
ii. Similarly, ii
~
~
C(ii)(u - ii),
u and hence u is the unique solution of
Equations of the form (5.7.1) arise as finite difference approximations to nonlinear partial differential equations as well as problems at resonance where the matrix is usually singular, irreducible, M-matrix. For such matrices the fol1owing result is known.
5.8.
151
Problems
Theorem 5.7.3. Let A be an n x n singular, irreducible, M-matrix. Then (i) A has rank (n - 1); (ii) there exists a vector u > 0 such that Au = 0; (iii) Av ~ 0 implies Av = 0; and (iv) for any nonnegative diagonal matrix D, (A + D) is an M-matrix and (A + D)-t is a nonsingular M-matrix, if d jj > 0 for some i, 1 ~ i ~ n. As a consequence of Theorem 5.7.3, if we suppose that A in (5.7.1) is n x n singular, irreducible, M-matrix, then (A + M) is a nonsingular M-
matrix. Therefore assumption (iv) of Theorem 5.7.1 holds since nonsingular M-matrix has the property that its inverse exists and is greater than or equal to zero. Furthermore, in some applications to partial differential equations, the function f in (5.7.1) is of special type, namely feu) = (ft(Ut),/zCU2),'" '/n(u n». We can find in this special case lower and upper solutions such that assumption (ii) of Theorem 5.7.1 holds whenever we have lim sup};(uJ Sig(uj) < 0 IuJ... co
for each i.
In fact, by Theorem 5.7.3, there exists a ~ > 0 such that Ker A = and hence we can choose a A > 0 so large that f(A~) ~
0
Span(~)
and f( -A~) ~ O.
Letting v=-A~, w=Ag, we see that v~w, Av~f(v) and Aw~f(w). Assumption (iii) is simply a one-sided Lypschitz condition and in assumption (5.7.5), B is a diagonal matrix.
5.8. Problems 5.1.
Show that p( G'(x*» = 0 where G(x) is a simple root of F(x) = O.
5.2.
Show that if F: R S
=x -
(F'(X»-l F(x) and x*
R has a Lypshitz derivative then
~
IIF(x) - F(y) - F'(x)(y -
x)11 ~ hllx - y112.
5.3.
Show that if in the Newton-Kantarovich theorem one supposes IIF'(x)-tll bounded for all x E R S , then [x, - x*II ~ cllxn - 1 - x*II2•
5.4.
Show that in the hypothesis of Theorem 5.3.4, the error satisfies the inequalities (5.3.12) and (5.3.13).
5.5.
Suppose that
x, E R
S ,
n
~ 0 and Il.Ilxnll ~
n
.Ilt
.Il
n- t
lI.1lxn-t1I, where
a positive converging sequence with to = 0 and lim tn n~CO
that
Xn
converges to a limit x".
= r:
tn is
Show
152
5.
Applications to Numerical Analysis
IlllxnI s; ( Iltn)'Y (1IIlxn-lll)Y Ilt
5.6.
As in the previous exercise suppose that
5.7.
Show that the solution of
5.8.
Find the constant k in the previous problem and deduce that for Newton's method one has:
Zo -
n- I
~z~
IS given by Zn = 1 - Zn 1 - (1 - 2zo) 1/ 2 cotght cz"), with k appropriately chosen.
'I
X
*-
Xn
Zn+1
I s; -(1 2 f3y
=
2'
2zo)1/2 - -(J2 - " 1 - (J
where 1 - Zo - (1 - 2Z0) 1/ 2 (J = ---=--'-----=..:--:-= 1 - Zo (1- 2Z0 ) 1/ 2 '
+
5.9.
Obtain the result of Theorem 5.3.3 by considering the first integral with Zo = ~. (Hint: in this case the first order equation becomes Zn+1 = ~(1 + zJ.)
5.10. Consider the iterative method
= JL[AG(x n ) + (1- A)G(xk, Xk-I)] + (1- JL)Xk, A, JL E [0,1] and G, G are defined as in Theorem 5.3.5. Show Xk+1
where that JL and A can be chosen such that the convergence becomes faster. 5.11. Solve equation (5.5.1) directly and show the growth of the errors. 5.12. Obtain the relations (5.5.2), (5.5.3) and (5.5.4). 5.13. Show that the sweep method is equivalent to the Gaussian elimination of the problem Ay = g when the problem 5.5.1 is stated in vector form, see (5.4.16).
5.9. Notes The discussions of Sections 5.1 to 5.3 have been introduced only to show examples of application of difference equations to numerical analysis. More details can be found in the excellent books by Ortega and Rheinboldt [127] and Ostrowsky [132]. Theorem 5.3.3 is a simplified version of the original one, see [129], [127], while the estimates given before 5.3.4 (b) seem to be
5.9.
Notes
153
new. For other estimates of the error bounds for the Newton's method and Newton-like methods, see also [41], [109], [111], [179], [144]. For the effect of errors in the iterative methods, see Urabe [168]. A very large source of material on the problem discussed in Section 5.4 can be found in Gautschi's review paper [52]. Theorem 5.4.1 has been taken from Olver [123]. A detailed analysis, in a more general setting, of the Miller's algorithm is given in Zahar [180]. Applications to numerical methods for ODE's can be found in Cash's book [23], while a large exposition of applications as well as theoretical results are in Wimp's book [56], see also Mattheji [104]. The "sweep" method can be found in the Godunov and Rayabenki's book [59], where it is also presented for differential equations. An improvement of the Olver's algorithm can be found in [169]. An application of the sweep method to the problem of Section 5.4 can be found in Trigiante-Sivasundaram [167]. The contents of Section 5.6 are adapted from [89, 126, 127], while the material of Section 5.7 is taken from [91], see also [85, 87, 89].
5.0.
Applications to Numerical Analysis
Introduction
Despite the fact that the two theories have been developed almost independently, the connections between numerical analysis and the theory of difference equations are several. In this chapter, we shall explore some of these connections. In Sections 5.1 to 5.3, iterative methods for solving nonlinear equations are discussed and the importance of employing the theory of difference inequalities is stressed. Sections 5.4 and 5.5 deal with certain classical algorithms from the point of view of the theory of difference equations. Sections 5.6 and 5.7 are devoted to the study of monotone iterative techniques, which offer monotone sequences that converge to multiple solutions of nonlinear equations. This study also includes extension of monotone iterative methods to nonlinear equations with singular linear part as well as applications to numerical analysis. In Section 5.8 we provide problems of interest.
5.1.
Iterative Methods
By iterative methods one usually means methods by which one is able to reach a root of a linear or nonlinear equation. Let us consider the problem of solving the algebraic equation of the form F(x) = 0,
(5.1.1)
127
128
5.
where F: R
S
~
R
S ,
Applications to Numerical Analysis
which is usually transformed into the iterative form:
Xno = Xo. (5.1.2) The function f(x) is called the iterative function and is defined such that the fixed points of f are solutions of (5.1.1). There are many ways to define the iterative function f. The essential criterion is that the root x* of (5.1.1), which we are interested in, is asymptotically stable for (5.1.2), although it may not be sufficient to ensure rapidity of convergence. This, however, will imply that if Xo lies in the region of asymptotic stability of x*, the sequence will converge. If f is linear, we know that asymptotic stability can be recognized by looking at the eigenvalues of A (the iterative matrix) and indeed this is what is done in the study of linear iterative methods such as Jacobi, Gauss-Seidel, SOR, etc. Moreover, in the linear case, asymptotic stability implies global asymptotic stability. For nonlinear equations the situation becomes more difficult. In this case there are essentially three kinds of results that one can discuss: Xn+l
= f(x n ) ,
A. LOCAL RESULTS. These results ensure asymptotic stability of x*, but they have nothing to say on the region of asymptotic stability. Usually theorems in this category start with the unpleasant expression "If Xo is sufficiently near to x* ...." B. SEMI LOCAL RESULTS. The results in this category verify that an auxiliary positive function (usually the norm of the difference of two consecutive terms of the sequence) is decreasing along the sequence itself. The sequence is supposed to (or one proves that it) lie in a closed set D c R Then one can infer that the auxiliary function has a minimum in D and this minimum is located at x*. We shall see that this usually requires that x* is exponentially stable. Connected to these types of results are those requiring the stronger condition of contractivity off in a closed set D c R Contractivity implies exponential stability of the fixed point x*. C. GLOBAL RESULTS. These results say that the sequence x, given by (5.1.2) is almost (except, on a set of measure zero) globally convergent. S
•
S
•
Of course, results in the class (C) are very few in the nonlinear case. One such result is due to Barna and Smale, which says that for a polynomial with only real root, Newton's method is almost globally convergent, which we shall not present here. We shall, however, discuss in the next two sections some of the most important results in the classes (A) and (B).
5.2. Local Results
"
Let us begin with the following main result. Here f: R x* = f(x*).
S
~
R and x* satisfies S
5.2.
129
Local Results
Theorem 5.2.1. Suppose that f is differentiable with Lipshitz derivative, and the spectral radius off'(x*) is less than one. Then x* is asymptotically stable. Proof. By the transformation eTj difference equation en + 1
= f(x n) - f(x*) = f'(x*)en + =
where
f'(x*)en
=
xTj - x*, from (5.1.2) we obtain the
L
{f'(x*
+ s(xn -
x") - f'(x*)} ds en
(5.2.1)
+ g(en)
l~~ II~~:II)II = 0,
(see problem 5.2). The equation (5.2.1) is in the
form required by Corollary 4.7.2, which assures the exponential asymptotic stability of x*. Note that the exponential asymptotic stability does not imply in general that /lxl - x*11 < Ilxo - x"], unless (see section 4.2) aT] < 1, where the constant a depends on the norm used. The result established in Chapter 4 can be used to prove the next theorem on the convergence of nonstationary iterative methods defined by (5.2.2) where q : N;o x D c R ~ R , XO E D. It is obvious that nonstationary iterative methods are, in our terminology, nonautonomous difference equations. • S
Theorem 5.2.2.
S
Consider the difference equation
yno = X o ,
(5.2.3)
where f: D ~ R S , and f is locally Lipshitzian. Suppose that y* E D is an asymptotic stable solution for (5.2.3) and that the solutions y(n, no, x o) of (5.2.3) and x(n, no, xo) of (5.2.2) are in D for n E N;o' Assume further that
Ilq(n, x) - f(x)11 ~ Lnllx where L; ~ 0 for n ~ suitably chosen.
Proof.
y*ll,
(5.2.4)
Then the solution x(n, no, xo) ~ y*, when
00.
Xo
is
Rewrite (5.2.2) as Xn+I
= f(x n ) + R(n, x,)
(5.2.5)
and consider (5.2.5) as the perturbed equation 4felative to (5.2.3) with (5.2.6)
130
5.
Applications to Numerical Analysis
as the perturbation, which tends to zero as x, 4.11.1. •
~
y*. Then apply Corollary
Other results of local type can be obtained by using the comparison principle given in section 1.6. For example: Theorem 5.2.3.
Let f: D c R S ~ R S continuous, g: R+ ~ R+. Suppose that
(1) there exists x*
E
Int D such that for all x
E
B(x*, ll) < D, II
> 0,
(5.2.7) Ilf(x) - x*11 :5 [x - x*11 + g(llx - x*Ii); u + g(u) is nondecreasing with respect to u and g(O) = 0;
(2) G(u) = (3) the trivial solution of
Un + 1 =
G(u n ) ,
Uo = Ilxo - x*ll,
is asymptotically stable, with U o in the domain of asymptotic stability, and u, :5 Uo for all i. Then x* is an asymptotically stable fixed point for the equation
(5.2.8) Proof. Let Xo E B(x*, ll). Then by Theorem 1.6.1, we have [x, - x*11 :5 u; and, because Un is decreasing to zero, all x, E B(x*, 8) and the sequence will tend to zero. •
The next two theorems are applications of LaSalle's Theorem 4.8.6. Here we study the convergence of secant and Newton methods on the real line. These two methods are widely used in numerical analysis to find roots of equations and their generalizations are almost countless. On the real line, these two methods are defined by (5.2.9) and (5.2.10) respectively. Iff(zk) ¥- f(Zk+I), the difference equation (5.2.9) is well defined. Suppose that f is differentiable and that a is the simple root of f( z) = O. Using the transformation (5.2.11) and the mean value theorem f(a
+ e)
= f'(a)
+ g(a,
e)e 2 ,
(5.2.12)
S.2.
Local Results
131
we write the new difference equation corresponding to (5.2.9) having the fixed point at the origin in the form (5.2.13) where _ g(a, en+l)en+ l - g(a, en)e n M( a, en, en+ l ) ( . / a + en+!) - /(a + en)
(5.2.14)
If a is a simple root and g( a, e) is continuous and bounded, then Mt;«, ek, ek-l) is continuous and bounded if ek, ek-l are small enough. The second order difference equation (5.2.13) can be transformed into a first order system by introducing the variables x, = en and Yn = en+ l • Then the system is
= Yn, Yn+l = M(a, xn,Yn)xnYn' Lyapunov function V(x, y) = Ix[q + Iylq, X n+ l
By using the
q
2':
LlV(xn, Yn) = -(1 -IM(a, x n, Yn)Ynlq)lxnl
1, one obtains
q,
which is negative if W(x,y);;;: (l-IM(a,x,y)yj) > O. Let D(1) = {(x,y)I(!xl q + Iylq)l/q < 1)}. Since y = 0 implies W(x, y) = 1 > 0, there exists some 1) > 0 such that W(x, y) 2': 0 in D( 1). If the initial values Zo and ZI are such that (x o, Yo) and (Xl> Yl) E D( 1), then (x n , Yn) will remain in D (see Corollary 4.8.5) and will converge to the maximum invariant set contained in the set E = {(x, y) E 15(1) Ix = O}. The only invariant set with Xl = 0 is the origin, so we obtain (xn, Yn) ~ (0,0) for n ~ 00. Corollary 4.8.5 can be used to study the convergence of (5.2.10). Because there is no more difficulty involved we shall consider the case where (5.2.10) is defined in R S • The root a, which we are interested in, is simple. The change of variable similar to (5.2.11) allows us to write the method in the form where M l ( en) -_ [d/(a It can be
+ en)]-l and M ( en) -_ [d/(a + en)] . 2
ax shown that for I enII <
IIMI(en)[Mz(en)en
ax
1), -
/(a
+ en)] II =:; k(1)lle~ll.
132
S.
Then, taking Vee)
= IIell, one
Applications to Numerical Analysis
can obtain
AV(e):s -(1- k( 17)llell)llell
;:: - W(llell),
showing that A V :s 0 for k( 17) II en II :s 1. Using Corollary 4.8.5, we get a region of convergence D( 170) = {x Illx - a II < 170} with 170 = mine 17, 1/ k( 17». The parameter 17 can be chosen to maximize 170 and consequently the region of convergence.
5.3. Semilocal Results Results of semilocal type are obtained either by requiring contractivity or some of its generalizations on a compact D c R or by requiring the weaker condition that some auxiliary function decreases over a sequence. We shall consider some cases of the latter type. We begin discussing the simplest of such results in some detail to bring out clearly the arguments used. S
Theorem 5.3.1. Suppose that all x, obtained by (5.1.2) lie in a compact Do c R S and that
Ilxn+2 - xn+111:s allxn+l - xnll, with a < 1. Then the sequence converges to x* Proof.
E
(5.3.1)
Do.
Setting (5.3.2)
we have, because of (5.3.1), We can now apply Theorem 1.6.1 to obtain solution of Un + 1
Yn:S Un'
where
Un
is the (5.3.4)
= aU n
provided that Yo :S Uo. Now the equation (5.3.4) has the solution u; = 0 exponentially asymptotically stable, which means (see Theorem 4.1.1) that I;:l uj is bounded. In 00
fact from (5.3.4) we get Ij=o uj =
1 I-a
--Uo.
It then follows that
Uj,j = 1,2, ... , if we choose Uo = Ilxl IIxn+p
-
Xn
II :S IIxn+p
-
Ilxj + 1 -
Xjll
:S
- xoll. Moreover, for p ~ 1, we obtain x n+p - 1 + x n+p - 1 - X n- p-2 + ... + X n+1 - X nII
p-l
:S j=O L Ilxn+j+l -
p-l
xn+J
:S j=O L un+j •
133
5.3. Semllocal Results
Because the right-hand side can be made less than an arbitrary positive quantity when n is chosen properly, this inequality shows that x, is a Cauchy sequence and thus must converge to x* in Do- From this result we can also find the distance between Xo and x*. In fact, from n-l
n-l
j=O
j=O
[x, - xoll ~ I IIXj+l - xjll ~ I
Uj>
it follows that
Ilx* -
xoll
~
1
00
I Uj j~O
=
--Uo· 1- a
(5.3.5)
Using this property one can avoid the a priori assumption that all x, lie in Do. Also, if we assume thatfis continuous, thenf(x*) = x*. Furthermore, if we assume more than (5.3.1), namely, x, y ED, Ilf(x) - f(y) II ~ k < 1, then the solution x* is unique in Do. • Using similar arguments, one can also prove the following result. Theorem 5.3.2. Suppose that xo, Xl lies in Do and Bt x«, L;'o Uj) c Do. Then the sequence converges to x* E B(xo,LJ:=o Uj)' The next theorems are on the convergence of Newton's method defined by (5.3.6) where F is the same function as in (5.1.1). The main result was first proved by Kantorovich and is usually called Newton-Kantorovich theorem. Theorem 5.3.3. ofR s,
Assume that, for X o, x, Y E D, where D is a convex subset
(l) F is differentiable in D; (2) IIF'(x) - F'(y)11 ~
(3)
II F'(xo)-lll
(4) f3'YTJ (5)
'Yllx - yll; 11F'(xo)-1 F(xo) II
~ f3,
=
TJ;
B( Xo, f31") c
D.
Then the sequence defined by (5.3.6) is well defined and converges to a root x* of F(x) = 0 contained in B(xo, I/f3'Y)· Proof.
For every
X E
Bt x«,
II {3'Y) we have
IIF'(x) - F'(xo)11 ~
'Yllx - xoll
< {3-I.
134
5.
Applications to Numerical Analysis
A bound of IIF'(x)-111 can then be obtained by using the Banach lemma (see Corollary A5.2), which yields <: IIF'(xo)-111 <: f3 x - I - IIF'(xo)-IIIIIF'(x) - F'(xo)ll- 1 - f3'Yllx - xoll'
IIF'( )-111
From this estimate and from hypothesis (4), it follows that for x E B(xo,I/f3'Y), F'(x) is nonsingular and the Newton's iterative function is defined. Moreover,
and (see problem 5.2) IIF(xn+1)11 = IIF(xn+l) - F(xn) - F'(xn)(xn+ 1 :s; hllxn+1 - x n11
2
-
xn)11
(5.3.8)
•
Substituting for F(x n + l ) in (5.3.7), we then find
IX
2
n+2 -
f3'Y Ilxn+1 - xnl1 Xn+1I :s; - 2 1 - f3'Y II.Xn+1 - Xo II'
Suppose for a moment that the sequence X n lies in B(xo, 1/ f3y), then we . gUn, ( "n "n-I) = f3'Yu~n ' The can apply Theorem 4.5.4 with L.j=1 Uj, L.j=1 U j
2(1 - f3y
Lu
j )
j=O
comparison equation is therefore f3Yu~
(5.3.9)
whose solution is Un = (,8'Y)-lr n- l. The origin is thus exponentially stable for (5.3.9). All x, lie in B(xo, 1/ ,81') because [x, - xoll < (f3'Y)-I2:;::11 2-
j
-
1
<
;1"
We may apply the
same arguments used
b ef ore. In f act, II xn+ p _ Xn II <: - "p-I L.j~O II xn+j+ 1 _ xn+j
II
<: -
"p-I _ (f3'Y)-1 L.j=O un+j 2
x 2:;::~ r n - j = (f3'Y)-lr n(l - r p ) , which shows that X n ~ x* and also gives the error estimate IIx* - xnll:s; (f3'Y)-lr n. The proof that x* is a root of F(x) follows from F(xn) = F'(Xn)(Xn+1 - x n) = [F'(xo) - (F'(x n) - F'(xO))](Xn+1 - x n) and then IIF(xn)ll:s; [IIF'(xo)11 + 'Yllxn - xoll]llxn+1 - xnll:s; [IIF'(xo)11 + 'Y(f3'Y)-I]llxn+1 - xnll. Taking the limit as n ~ 00 one has, by the continuity of F, F(x*) = O.
5.3.
Semilocal Results
135
Hence by Theorem 4.5.4, the conclusion follows. • The conclusion of Theorem 5.3.3 has been easily obtained because the solution of (5.3.9) can be given in the closed form. More refined estimates can be obtained if we assume in (5.3.9) that Uo < !(f~y)-l. In this case, the solution of (5.3.9) is more difficult to compute but it can be done. In fact, let us define the new sequence Zn by dZn- 1 = {3YUn, L l = O. Then . (5.. 3 9) IS . transtorme .. d'mto Zn+l = Zn + -1 (zn - Zn_l)2 , W h'IC h equation 2
1 - z;
possesses the first integral, that is, the sequence {z.} satisfies the first order . Zo - !z~ . of equation Zn+l = , where Zo = {3yuo . (See (1.5.12).) The solution 1 - Zn this equation is z; = 1 - (1 - 2Z0 ) 1/ 2 coth k;", (see problem 1.27) where k is a constant depending on the initial condition. (See problem 5.8 for its determination and for a more refined bound.) This solution exists if Zo s ! and then Uo s !({3y)-l. The case of equality has already been examined. Suppose that Ul < !({3y)-I. Then it follows that lim z; = 1 - (1 - 2Z0 ) 1/ 2 == z* and since for n~OO
• , we obtain z; - z* = 2(1 - 2z o) 1/2 e" k2 + . • Zn - z* Consequently, we get !~~ (Zn-l _ z*? = 1, which shows that the order of
large n, coth k2 n
= 1 + 2 e" k2
n+l
n
1
convergence of the sequence z.,, to the limit point, is quadratic. From
the
comparison
_l_(Zn - Zn-l) and hence {3y
[x, of z.;
xoll s l:;':-~ Uj S
Ilx
n -
principle
it
x*11 < ~(1 {3y
follows
that
Ilxn+1 -
Xn
II <
1
2Z0 ) 1/ 2 e- k 2 n + • Also for n > 0,
_l_ Zn_ 1 s·_l- z * because of the increasing behavior
{3y
(3y
•
Summing up these considerations we see that we have proved a variant of Theorem 5.3.3, which we state below. Theorem 5.3.3(a).
Assume that the hypotheses (1), (2), (3) of Theorem 5.3.3
hold. Suppose further that (4) f3yY]
Z*)
and (5) B( Xo, f31 c: Dare y satisfied. Then the sequence x, is well defined and converges quadratically to
x* contained in B( xo, (31y z*). Newton's method is invariant under affine transformation, that is, it is invariant under the transformation G = AF where A is an invertible matrix,
136
S. Applications to Numerical Analysis
whereas the hypotheses in the Newton-Kantarovich theorem are not. In fact, one has 0- 10 = F-1A-1AF = F- 1 F, while for example, the hypothesis (3) is not invariant under the same transformation, The following variant of Theorem 5.3.3 takes care of this situation. Theorem 5.3.3(b). and suppose that
Assume that the hypotheses (1) of Theorem 5.3.3 holds
II F'(XO)-I F(xo) II :s; a; IIF'(xo)-I(F'(y) - F'(x»11 :s; wllx - yll;
(2) (3) (4) (5)
aw:s;!; B(xo, w-1z*)
E
D, where z* = 1 - (1- 2wa)I/Z:s; 1.
Then the sequence x, remains in B(xo, w-1z*) and converges quadratically to a root x* of F(x).
Proof.
We have
Xn+Z - Xn+1 = -F'(Xn+I)-IF(xn+l)
= -F'(Xn+I)-I[F(xn) + F(xn+ 1 )
-
F(xn)]
-F'(Xn+I)-I{ F(x n) + F'(Xn)(Xn+1 - x n)
=
+ = -
f
[F'(x n + t(xn+1 - x n» - F'(x n)] dt(x n+1 - x n)}
F'(X n+1)-1 F'(xo){I F'(XO)-I[F'(xn + t(xn+1
-
»
xn
- F'(xn)] dt(x n+1 - Xn)}, from which it follows that Ilxn+z - xn+lll -s IIF'(xn+I)-1 F'(xo) II . ~ Ilxn+1 - x,
liZ,
Since F'(XO)-IF'(xn+ l ) == F'(XO)-I F'(Xn+l) - F'(xo)-I F'(xo) + F'(XO)-I F(xo), using hypothesis (3), we get 11F'(xO)-I(F'(xn+l) - F'(xo» II < wllxnn - xoll < z*:s; 1, if Xn+l E B(xo , w-1z*). Then, by the Banach lemma it follows that F'(x n+ l) is invertible and II F'(Xn+I)-1 F'(xo) II :s; 1-
II 1
W Xn + 1 -
II' Finally, we arrive at
Xo
IIXn+Z -
w Ilxn+1 - xnf xn+lll :s; -2 1 _ I _ II' w Xn+1 Xo
137
5.3. Semilocal Results
which is the same as before with f3y replaced by w. All the previous results then apply, and we are done. • The next result also shows the familiar quadratic decay of the errors in the Newton method, but it requires more assumptions on the inverse ofthe Jacobian.
Assume that F: Do c R S Moreover, suppose that
Theorem 5.3.4. (1) (2) (3) (4) (5)
~
R S and Do c R S is a convex set.
for x E Do, F is differentiable, for x E Do, F'(x) is invertible and IIF'(x)-111 for x, y E Do, 11F'(x) - F'(y)11 :5 yllx - yll, for Xo E Do, IIF'(xo)-IF(xo)II:5 TJ, and (\' = 4f3YTJ < 1.
:5
f3,
Then the Newton iterates remain in ei«; ro), where ro converge to a solution x* of F(x) = O. Proof.
= TJ 2:~~0 (\'ZH
and
From
IIx n+z -
xn+III =
IIF'(xn+ I)-IF(x n+I)11
:5
f3IIF(xn +1 ) II
and from the estimate of II F(xn + l ) I obtained in the previous theorem, we get (5.3.10) The comparison equation is therefore
f3y
Z
Un + 1 = TUn,
Uo
=
(5.3.11)
TJ,
whose solution is u; = TJ(\'zn_ l • Now applying comparison Theorem 1.6.1 we get (see problem 5.4)
[x, - xn+11 1:5
Un
= "1(\'zn_ l.
(5.3.12)
Since (\' < 1, the trivial solution 01 (5.3.11) is exponentially stable and hence it follows that x, is a Cauchy sequence. Furthermore, we also have [x, xoll :52:';:0 Uj = ro, which means that all the iterates remain in Bt x«, ro) and converge to x*. • Corollary 5.1.1.
The error [x" - Xn II satisfies the inequality
Ilx* - x, I :5 en IIxn -
Xn-I
liz,
for n
= 1,2, ....
(5.3.13)
As an application of Theorem 1.6.7, let us consider the two-step iterative method
138
5.
Applications to Numerical Analysis
Theorem 5.3.5. Suppose that G: D x D c R x R set Do c R the following inequality holds true S
S
~
R S and on some closed
S
,
IIG(x,y) - G(y,
z)ll:5 allx - yll + ~lly - zll, let there exist x o, XI E D such
where a + ~ < 1. Furthermore, that Xk E Do for k e: 1. Then the iteratives converge to the unique fixed point of G(x) = G(x, x). Proof. Setting Yk = Il xk + l - xkll, we obtain Yk+l:5 aYk + ~Yk-l' By Theorem 1.6.7, we then get Yk :5 Uk, where Uk is the solution of Uk+l = auc + ~Uk-I' which is Uk = coz~ + CIZ~, where ZI and Z2 are the roots of Z2 - az - ~ = O. Because of the assumption on a and ~, the two roots ZI and Z2 are less than 1 in modulus, which means that the zero solution of the comparison equation is exponentially stable. Using similar arguments as in the previous theorems, we can conclude that the sequence Xk is a Cauchy sequence. Since
IIG(x*) - Xk+111 :5IIG(x*, x") - G(x*, Xk-I)II + IIG(x*, Xk-I) - G(Xk' Xk-I)II :5 ~llx* - xk-t11 + allx* - xkll, the convergence of Xk shows that Xk+l ~ G(x*). The uniqueness of x* now • follows very easily and the proof is complete. The effect of perturbation on the iterative methods are of two different kinds. If the perturbations are small and tend to zero as n tends to 00, then the theorems on total stability will ensure that the fixed points continue to have asymptotic stability. In numerical analysis, however, perturbations may remain bounded for all n (roundoff errors, for example). In this case, a more convenient concept is the practical stability. If an iterative procedure is practically stable, the sequence of iterates will not tend to the solution x* but to a ball B(x*, B) surrounding x*. Inside B, the solution may oscillate but what is important is that it never leaves B. Of course it would be nice to have B as small as possible. Let us consider for example the iterative method
= Xo'
(5.3.14)
Ilf(x) - f(y) II :5 allx - yll
(5.3.15)
Xn+ 1 = f(x n ) ,
We assume that for x, y
E
Do c R
x na
S ,
with a < 1. Suppose that the errors in the computations perturb (5.3.14) by a bounded perturbation R n , that is
xn + 1 = f(x n ) + R n ,
xna = xo,
(5.3.16)
5.4.
139
Unstable Problems
with X n , xn E Do for all n. Then Theorem 4.11.2 can be applied. The conclusion is that the difference of the two solutions [x, - xn II will not exceed Rj(l- a), where R
B(Xn,~) I-a
=
sup Rn- The condition that all the balls
must be contained in Do may be very restrictive. In fact if a
is close to one, it follows that these balls are very large and the method is considered a bad one.
5.4.
Unstable Problems: Miller's, Olver's and Clenshaw's Algorithms
The situation that we shall analyze in this section and in the next is in a certain sense the opposite of that treated in the previous ones. In fact, there the limit point was important (and indeed in the theory of iterative processes the "solution" is the limit point), where the intermediate values z; are considered an unavoidable noise. Here the limit point will not be important (generally will not exist) and the important part ("the solution") will be a finite subset of the sequence itself. These problems are very typical in that part of numerical analysis which concerns the study of special functions, orthogonal polynomials, quadrature formulas, and numerical methods for ordinary differential equations. Consider the homogeneous scalar linear equation (see (2.1.5)) of order k. (5.4.1) It has k linearly independent solutions/l(n)'/in), ... '/k(n). Suppose that the solution Yn = 0 (n ~ 0) is unstable, and
. II(n) lim 1"( n ) = 0,
i = 2,3, ... , k.
(5.4.2)
n~oo.li
When this happens, the solution II(n) is said to be minimal. The problem is how to find II or a multiple of it. Even if we know the exact initial condition that generates the minimal solution, small errors in the calculations will, as usual, lead us to solve a perturbed equation LYn = en, whose solution (see (2.1.11)) will contain all the j;(n) and this will destroy the process (see for example [52], [175]), because the h(n) (i ~ 2) grow faster than II(n). The same problem appears, of course, in the nonhomogeneous case, where it is not excluded a priori that there exists a bounded solution even if all the hen) are unbounded. One is often interested in following this solution.
140
5.
Applications to Numerical Analysis
Consider, for example, the following first order equation Yn+l - (n
+ 1)Yn
= -1,
Yo = e.
(5.4.3)
The solution of the homogeneous part is Yn = en!, while the complete solution is Yn = n!( e -
2:.;=oh) ~ 2:.:=1 (n +l s
) (S) ,
which is bounded for
n ~ 00. If one tries to follow this solution, a small error in the initial data e (always existing because it has infinite digits), will be amplified by the factor n!. A simple idea is to find the way to look for the solution of the problem backwards (for n ~ -00), because in this case the origin becomes asymptotically stable for the homogeneous equation and the errors will be damped. This idea indeed works, as we will show in a simple case, which is the original proble~ to which it was applied (see [52] for references). Consider (5.4.6) where an, en ¥- 0 for all n, and the nonhomogeneous problem LYn = gn'
(5.4.7)
Suppose that two linearly independent solutions ¢n and l/Jn of the homogeneous equation are such that ¢n is minimal, that is
· ¢n 0. IIm-= "_00 t/Jn
(5.4.8)
Of course all solutions that are multiple of ¢n are minimal and vice versa. That means that only one appropriate initial condition is needed to determine the minimal solution we are interested in. The general solution of the nonhomogeneous problem is then (5.4.9) where
Yn
is a particular solution that will be supposed minimal too; that is,
· Yn 0 IIm-= I/Jn '
n_OO
(5.4.10)
and Yo = O. We are interested in the solution Yn
=
Yo,/, +_ ¢o't'n Yn'
A class of methods are based on the following result.
(5.4.11)
5.4.
141
Unstable Problems
Theorem 5.4.1.
Consider the boundary value problem Ly
= gn,
y~N)
= Yo,
Y<:')
= 0,
(5.4.12)
where N is a positive integer, and suppose that conditions (5.4.8) and (5.4.10) are verified. Then the problem (5.4.12) has a solution and moreover, for fixed n, y~N) ~ Yn as N ~ 00.
Proof. Since Yn is a particular solution of the nonhomogeneous equation, any other solution of the same equation will be y~N)
= c~N)cPn + ciN)ljJn + Yn'
(5.4.13)
The boundary conditions give for C~N) and ciN) the values YN C(N) _ I
-
IjJN
ljJo
It follows that C~N) and ci
(N) __ cPo (N)
.
cPo _ cPN'
C2
-
ljJo
CI'
IjJN
N)
tend to zero as N ~ 00 and the claim follows. Note that the condition Y<:') = 0 can be replaced by Y<:') = k where k is any fixed value (for example an approximation of YN if available). In the homogeneous case y~N) can be obtained starting from j<:') = 0, j~NlI = 1, and then using the difference equation backwards ~(N)
Yn-I
an _
b., :
c,
en
= --Yn+1 --Yn
obtaining a value y~N) different from Yo. Multiplying then the sequence by the scale factor Yolj~N), one obtains y~N). This is the Miller's algorithm. A disadvantage of this algorithm is that one does not know a priori which value of N must be used. A modified version of it, which works in the nonhomogeneous case is also obtained by considering that the boundary value problem (5.4.9) is equivalent to the system of equations Ay(N) = b, (5.4.14) where
bl
0
0
al
CI
A=
0
0 aN-I
0
0
bN - I
CN-I
and Y
(N) _
-
b
«(N)
Yt
= (gt
(N)
"'YN-t
)T
,
- cIYo, ... , gN-t) T.
142
5.
Applications to Numerical Analysis
Some clever method to solve this system can be used to avoid the growth of the errors (see [23], [24]). By solving the system in an appropriate way, it is also possible to determine the best value of N (Olver's algorithm). • Similar in motivation to the previous algorithms and indeed related to them is the following, known as Clenshaw's algorithm. Here the problem is to evaluate the sum L~=o bnun , where b; is a given sequence and u.; is supposed to satisfy a k th linear difference equation LUn = O. The algorithm uses the result stated in Theorem 2.1.8, that is, one considers the transposed equation LTYn = bn, YN+j = 0 (j = 1,2, ... , k), obtaining N
I
n=O
k-l
b.u;
=I
j=O
j
Yj
I
i=O
Pi(j - k)Uj_i'
Suppose, for example, we want to compute IN = Ln~o bnTn(z), where Tn(z) are Chebyshev polynomials satisfying the equation Tn + 1 - Zz'T; + T; = O. The transpose equation is N
Yn - 2ZYn+l YN+1
+ Yn+2
=
bn,
n
= N, N
- 1, ... , 0
= YN+2 = 0,
which can be solved recursively obtaining Yl and Yo. From the quoted result it then follows IN = Yo + Yl [T1(z) - 2zTo(z)] = Yo - ZYI' It is interesting to note that the sum can be obtained without the knowledge of the polynomials Tn(z) except for To(z) and T1(z) and also the reduction of the operations involved.
5.5.
Unstable Problems: the "SWEEP" Method
The topic discussed in this section is very similar to the previous one, but historically it has different sources. It arises typically in that part of numerical analysis that concerns the approximation of solutions of ODE and PDE. We shall consider, for simplicity, the autonomous case. The more general case is not, in principle, more complicated. Let us consider Yn+l - 25Yn + Yn-l
= gn,
y(O) = a,
Yk
= {3,
(5.5.1)
with 5 > 1, k> O. This problem can be solved in many different ways. It can be solved rewriting it as a system of k - 1 equations similar to (5.4.14). In order to control the growth of errors it is better to treat the problem as a difference equation.
5.5.
143
Unstable Problems
The characteristic polynomial of (5.5.1) has two positive roots, one less than one and the other greater than one. This means that the origin is not asymptotically stable for the homogeneous part and this, as usual, gives troubles because the errors are amplified. One way to avoid this difficulty is to transform the linear second order equation (5.5.1) into nonlinear first order equations (see problem 5.12). 1 Xn+1 = 2i> _ x ' (5.5.2) n
and (5.5.3) Furthermore, one verifies that Yn-l
= XnYn + z.:
(5.5.4)
The solution is then obtained by computing recursively the sequences {x n } , {zJ using (5.5.2) and (5.5.3) with XI = 0 and ZI = a up to n = k - 1 and then computing {Yn} using (5.5.4) from n = k to n = 1 with Yk = (3. The advantage is that (5.5.2) has a limit point asymptotically stable and the sequence x, converges monotonically to this point, while the sequences (5.5.3) and (5.5.4) remain bounded.
Theorem 5.5.1. If Xo = 0, 8> 1, the sequence {x n } converges monotonically to the first root p of the characteristic polynomial of (5.5.1) and, moreover, 0:$ x, < p. The critical points of (5.5.2) are the roots of the equation x 2 - 2i>x + 1 = 0, which is the characteristic polynomial of (5.5.1). Let p be the root less than one. Changing the variable h; = p(p - x n ) , equation (5.5.2) becomes Proof.
p2 hn + 1 = --h-hn • 1+ n
(5.5.5)
One shows at once (for example using the Lyapunov function Vn = h~) that the half line h > 0 is contained in the asymptotic stability region of the origin, and moreover 0 < h; :$ p2, from which it follows that -p :$ x, - P < O.
•
(5.5.6)
The next theorem is not really necessary in this case, but it could be useful in the more general case.
s.
144
Applications to Numerical Analysis
Theorem 5.5.2. In the hypothesis of Theorem 5.5.1 the solution x, = P is practically stable for (5.5.2). Proof.
The conditions of Theorem 4.11.2 are satisfied for x, < p.
•
The last result ensures that for equation (5.5.2) the errors will not be amplified. Let us consider the equations (5.5.3) and (5.5.4). Theorem 5.5.3. In the hypothesis of Theorem 5.5.1, the sequence {z.} remains bounded, even if the two equations are perturbed with a bounded perturbation. Proof.
From
p2 - 20p
+ 1 = 0, it follows that 20' - p = 1/ p and = P - gn( + Zn ) and f rom (5.5.6 )
. ( ) can be wntten 5.5.3 as Zn+1
IZn+tl:::; p(lgnl + IZnl)·
I - p xn - p
Using Corollary 1.6.1, it follows that IZnl:::; pna + L;:~ pn-s+1lgsl from which we see that if'[g.] is bounded z; is bounded. Clearly if we add to gn a bounded perturbation Tn (any sort of errors), the statement remains true. • Similar arguments can be applied to (5.5.4) since it follows that IYn-tl :::; plYnl + IZnl· The method just described is known as the "SWEEP" method (see [28]).
5.6. Monotone Iterative Methods In this section, we shall develop a general theory for a broad class of monotone iterations. This class of iterations includes Newton's method as well as a family of methods, which are called Newton-Gauss-Seidel processes that are obtained by using the Gauss-Seidel iteration on the linear systems of Newton's method. As before, we are interested in finding the solutions of
0= f(u), where f
E
C[R
n
;
R
n
].
(5.6.1)
Let us first split the system (5.6.1) as
0= j;(u i , [u]p" [u]q,),
(5.6.2)
where, for each i, 1:::; i:::; n, Pi + qj = n - 1 and u = (Ui, [u]p" [u]q')' Let v, W E R" be such that v :::; w. Then v, ware said to be coupled quasi lower and upper solutions of (5.6.1) if
o:::;j;(v [v]p" [w]q,), o-z.j;(w j, [w]p" [v]qJ i,
145
5.6. Monotone Iterative Methods
Coupled quasi extremal solutions of (5.6.1) can be defined easily. Vectorial inequalities mean the same inequalities between the components of the vectors. We are now in a position to prove the following result.
Theorem 5.6.1. Assume that f E C[R n , R n ] and possesses a mixed quasimonotone property, that is, for each i,};(uj, [u]p" [u]q') is non decreasing in [u]p, and nonincreasing in [u]q,. Suppose further v, ware coupled quasi lower and upper solutions of (5.6.1) and
j;(Uj, [u]p" [u]q) - };(iij, [u]p" [u]q.) 2:= -Mj(uj - iij) whenever v ~ ii -s U ~ wand M, > O. Then there exist monotone sequences {v n }, {w n } such that Vn ~ p, Wn ~ r as n ~ ex) and p, r are coupled minimal and maximal solutions of (5.6.1) such that v ~ p ~ u ~ r ~ wfor any solution u of (5.6.1). Proof.
For any
771,772 E [v, w] = [u E R" : v ~ u ~ w].
Consider the
system
u, = Mi l};( 77ti, [77I]P,' [772]q.)
+ 77ti,
i
= 1,2, ... , n.
Clearly U can be uniquely defined, given 771, 772 E [v, w]. Therefore, we can define a mapping A such that A[ 771, 772] = U. It is easy to show that A satisfies the properties: (i) 0 ~ A[v, w], 02:= A[w, v]; (ii) A is mixed monotone on [v, w]. Then the sequences {v n }, {w n } with Vo
= v, Wo = w can be defined as follows
Furthermore, it is clear that {v n } , {w n } are monotone sequences such that v ~ V n ~ W n ~ wand consequently lim V n = p, lim W n = r exist and satisfy the relations n_OO n_OO
By induction, it is also easy to show that if (Ul> U2) is a coupled quasi solution of (5.6.1) such that v ~ Ul> U2 ~ w, then V n ~ Ul> U2 ~ W n and consequently (p, r) are coupled quasi extremal solutions. Since any solution U of (5.6.1) is a coupled quasi solution, the conclusion of the theorem follows and the proof is complete. • If f does not possess the mixed quasi monotone property, we need a different set of assumptions to generate monotone sequences that converge to extremal solutions. This is the content of the following results.
146
S.
Theorem 5.6.2.
Applications to Numerical Analysis
Assume that
(i) there exist v, W E R" with v :s; W such that 0 :s; f( v), 0 =:: f( w); (ii) there is a n x n matrix M such that f(y) - f(x) =:: -M(y)(y - x), whenever v :s; x :s; y :s; w. Then the sequence {w n}, given by Wn + 1
= Wn + Bnf(wn ) ,
(5.6.3)
where B; is any nonnegative subinverse of M( wn), is well defined, and {w n} is monotone non increasing such that lim Wn = r. If there exists a nonsingular n~OO
matrix B =:: 0 such that lim inf B; =:: B, then r is a maximal solution of (5.6.1). n-e co
Set v = Vo and W = Wo' From B o > 0 and f( wo) :s; 0 it follows that Using the fact that B o is a subinverse of M (w o), we find for any [v, w],
Proof.
WI :s; W o'
u
E
u
+ Bof(u) =
WI -
:s; WI -
Hence, in particular, vo:S; Vo f( WI)
:s;
(w o - u) - Bo[f(w o) - f(u)J
[I - BoM(wo)](wo - u):s;
+ Bof( vo) :s;
f( wo) + M( woH Wo -
WI) =
WI'
WI'
Similarly, we obtain
[I - M( wo)BoJf(wo) s
o.
Proceeding similarly we see by induction that f(w n ) =:: 0,
n = 1,2, ....
(5.6.4)
Consequently, as a monotone nonincreasing sequence that is bounded below, {w n } has a limit r 2: Vo. If u is any solution of (5.6.1) such that u E [v, w J, then u = u + Bof( u) :s; WI' then by induction u :s; W n for all n. Hence u :s; r. Finally, the continuity of f and the fact lim inf B; =:: B, where B =:: 0 is nonsingular, (5.6.4) yields n
= -(lim inf B; )f( r) =:: -
Bf( r)
2:
0,
n~OO
which implies f(r) = 0 completing the proof.
•
By following the argument of Theorem 5.6.2, one can prove the next corollary.
5.7.
147
Monotone Iterative Methods
Corollary 5.6.1.
Let the assumptions of Theorem 5.6.2 hold. Suppose that M(y) is monotone non increasing in y. Then the sequence {v n } with Vo = v, given by V n+1
= Vn + Bnf(vn)
is well defined, monotone nondecreasing such that lim minimal solution of (5.6.1). n-+OO
Vn
= p, and is the
The case of most interest is when p = r, because then the sequences {v n }, {w n } constitute lower and upper bounds for the unique solution of (5.6.1). The following uniqueness result is of interest.
Theorem 5.6.3. Suppose that f(y) - f(x):5 N(x)(y - x), where N(x) is nonsingular matrix and N(X)-l ~ O. Then if (5.6.1) has either maximal or minimal solution in [v, w], then there are no other solutions in
[v, w]. Proof. Suppose r E [v, w] is the maximal solution of (5.6.1) and u is any other solution of (5.6.1). Then
0= f(r) - f(u)
:5
E [v, r]
N(u)(r - u).
Since N(U)-l ~ 0, it follows that r ~ u and hence r = u. A similar proof holds if the minimal solution exists. The proof is complete. •
5.7.
Monotone Iterative Methods (Continued)
Consider the problem of finding solutions of Ax = f(x),
(5.7.1)
where A is an n x n matrix and f E C[R n, R n ] , which arises as finite difference approximation to nonlinear differential equations. If A is nonsingular, writing F(x) = f(x) - Ax = 0, one can study existence of multiple solutions by employing the method of upper and lower solutions and the monotone iterative technique, described in Section 5.6. In this section we extend such existence results to equation (5.7.1) when A is singular. For convenience let us split the system (5.7.1) and write in the form (5.7.2)
148
5.
Applications to Numerical Analysis
where for each i, 1 ~ i ~ n, u = (Ui' [U]Pi' [u]q') with Pi + qj = n - 1 and (Au); represents the i l h component of the vector Au. As before, a function f E C[R", R"] is said to be mixed quasi-monotone if, for each i, /; is monotone non decreasing relative to [u ]p, components and monotone nonincreasing with respect to [u ]q, components. Let v, w E R" be such that v ~ w. Then v, ware said to be coupled quasi lower and upper solutions of (5.7.1) if
Coupled quasi-extremal solutions and solutions can be defined with equality holding. We are now in a position to prove the following result. Theorem 5.7.1. Assume that (i) f E C[R n , R n ] and f is mixed quasi-monotone; (ii) v, ware coupled quasi lower and upper solutions of (5.7.1); (iii) /;(u;, [u]p" [u]q') - /;(u;, [u]p" [u]q,) ~ -Mj(u j - uJ whenever v ~ u ~ u ~ wand M, > 0, for each i; (iv) A is an n x n singular matrix such that A + M = C is nonsingular where M is a diagonal matrix with M, > 0 and C- I ~ O. Then there exist monotone sequences {v n }, {w n } such that V n ~ p, W n ~ r as n ~ 00 and p, r are coupled quasi-extremal solutions of (1.1) such that if u is any solution of (1.1), then v :5 P ~ u ~ r ~ w.
Proof. For any YJ, I-t system
E
[v, w]
=
[x
E
Rn:v
~
x
~
w], consider the linear
(5.7.3)
Cu = F( YJ, I-t),
where C=A+M and for each i, Fj(YJ,I-t) =/;(YJk, [YJ]p,,[I-t]q,) + MjYJj. Clearly u can be uniquely defined given YJ, I-t E [v, w] since the matrix C is nonsingular. Furthermore, we have Av ~ F(v, w), Aw ~ F(w, v), and F is mixed monotone. Consequently we can define a mapping T such that T[ YJ, I-t] = u and show easily, using the fact that C- I ~ 0, that (a) v ~ T[v, w], w ~ T[w, v]; (b) T is mixed monotone on [v, w]. Then the sequences {v n } , {wn } with follows:
Vo
= v,
Wo
=w
can be defined as
S.7. Monotone Iterative Methods
149
Furthermore, it is evident from the properties (a) and (b) that {vn } , {w n } are monotone such that Vo ~ VI ~ ••• ~ vn ~ Wn ~ ••• ~ WI ~ Wo0 Consequently, lim e, = p, lim Wn = r exist and satisfy the relations n-JoOO
n-JoCO
By induction, it is also easy to prove that if (u l , U2) is a coupled quasisolution of (5.7.1) such that v ~ UI, U2 ~ w, then V n ~ ul , U2 ~ Wn for all nand hence (p, r) are coupled quasi-extremal solutions of (5.7.1). Since any solution u of (5.7.1) is a coupled quasi-solution, the conclusion of the theorem follows and the proof is complete. • The case of most interest is when p
=
r, because then the sequences {vn },
{w n } constitute lower and upper bounds for the unique solution of (5.7.1).
The following uniqueness result is thus of interest. Corollary 5.7.1. for each i, that
If, in addition to the hypotheses of Theorem 5.7.1, we assume (5.7.5)
where B is an n x n matrix such that (A - B) is nonsingularand (A - B)-I O. Then u = p = r is the unique solution of (5.7.1) such that v ~ u ~ w. Proof.
Since p
~
r, it is enough to prove that p
2::
2::
r. We have by (5.7.5)
[A(r - p)]j ~j;(ri' [r]pi' [p]q) - j;(Pj, [P]Pi' [r]q) ~ [B(r - p)]j.
Consequently, (A - B)(r - p)
~
0, which implies r ~ p.
•
Remark 5.7.1. If qi = 0 for each i, then Theorem 5.7.1 shows that r, pare maximal and minimal solutions of 5.7.1 and Corollary 5.7.1 gives the unique solution.
If f does not possess the mixed quasi monotone property, we need a different set of assumptions to generate monotone sequences that converge to a solution. This is the content of the next result. Theorem 5.7.2. Assume that (iv) of Theorem 5.7.1 holds. Further, suppose that V, W E R" such that v ~ W, Av
~f(v)
- B(w - v),
-B(x - y)
~f(x)
Aw 2::f(w) - f(y)
~
+ B(w -
B(x - y)
v),
(5.7.6) (5.7.7)
150
5.
Applications to Numerical Analysis
whenever v ~ y ~ x ~ w, B being an n x n matrix of nonnegative elements. Then there exist monotone sequences {v n }, {w n } that converge to a solution u of (5.7.1) such that
provided that (A - B) is a nonsingular matrix. Proof.
We define F(y, z)
= ![f(y) + fez) + B(y -
z)).
(5.7.8)
It is easy to see that F(y, z) is mixed monotone and
-B(z - z)
~
F(y, z) - F(y, z)
whenever z, z, y, Y E [v, w] and i
~
~
B(y - y)
(5.7.9)
z, Y ~ y. In particular, we have
F(y, z) - F(z, y)
= B(y -
z).
(5.7.10)
From (5.7.8) we obtain -B(w - v) ~ F(v, w) - f(v), which yields Av ~ f(v) - B(w - v) ~ F(v, w).Similarly Aw ~ F(w, v). Finally, it follows from (5.7.8) that F(x, x) = f(x). Consider now for any Tl, JL E [v, w], the linear system given by Cu = O( Tl, JL) where C = A + M, M being the diagonal matrix with M; > 0 and for each i, OJ( Tl, JL) = F; ( Tl, JL) + MjTli' Proceeding as in proof of Theorem 5.7.1, we arrive at Ap = F(p, r), Ar = F(r, p). Using (5.7.9), we see that A(r - p) = F(r, p) - F(p, r) = B(r - p) and this implies u = r = p is a solution of (5.7.1). The proof is complete. • Corollary 5.7.2. If in addition to the hypotheses of Theorem 5.7.2 we also have f(x) - fey) ~ C(y)(x - y) for x, y E [v, w], where C(y) is an n x n matrix such that [A - C(y)] is nonsingular and [A - C(y)r 1 ~ 0, then u is the unique solution of (5.7.1).
Proof.
If ii is another solution of 5.7.1, we get A(u - ii) = feu) - f(ii)
which yields u (5.7.1). •
~
ii. Similarly, ii
~
~
C(ii)(u - ii),
u and hence u is the unique solution of
Equations of the form (5.7.1) arise as finite difference approximations to nonlinear partial differential equations as well as problems at resonance where the matrix is usually singular, irreducible, M-matrix. For such matrices the fol1owing result is known.
5.8.
151
Problems
Theorem 5.7.3. Let A be an n x n singular, irreducible, M-matrix. Then (i) A has rank (n - 1); (ii) there exists a vector u > 0 such that Au = 0; (iii) Av ~ 0 implies Av = 0; and (iv) for any nonnegative diagonal matrix D, (A + D) is an M-matrix and (A + D)-t is a nonsingular M-matrix, if d jj > 0 for some i, 1 ~ i ~ n. As a consequence of Theorem 5.7.3, if we suppose that A in (5.7.1) is n x n singular, irreducible, M-matrix, then (A + M) is a nonsingular M-
matrix. Therefore assumption (iv) of Theorem 5.7.1 holds since nonsingular M-matrix has the property that its inverse exists and is greater than or equal to zero. Furthermore, in some applications to partial differential equations, the function f in (5.7.1) is of special type, namely feu) = (ft(Ut),/zCU2),'" '/n(u n». We can find in this special case lower and upper solutions such that assumption (ii) of Theorem 5.7.1 holds whenever we have lim sup};(uJ Sig(uj) < 0 IuJ... co
for each i.
In fact, by Theorem 5.7.3, there exists a ~ > 0 such that Ker A = and hence we can choose a A > 0 so large that f(A~) ~
0
Span(~)
and f( -A~) ~ O.
Letting v=-A~, w=Ag, we see that v~w, Av~f(v) and Aw~f(w). Assumption (iii) is simply a one-sided Lypschitz condition and in assumption (5.7.5), B is a diagonal matrix.
5.8. Problems 5.1.
Show that p( G'(x*» = 0 where G(x) is a simple root of F(x) = O.
5.2.
Show that if F: R S
=x -
(F'(X»-l F(x) and x*
R has a Lypshitz derivative then
~
IIF(x) - F(y) - F'(x)(y -
x)11 ~ hllx - y112.
5.3.
Show that if in the Newton-Kantarovich theorem one supposes IIF'(x)-tll bounded for all x E R S , then [x, - x*II ~ cllxn - 1 - x*II2•
5.4.
Show that in the hypothesis of Theorem 5.3.4, the error satisfies the inequalities (5.3.12) and (5.3.13).
5.5.
Suppose that
x, E R
S ,
n
~ 0 and Il.Ilxnll ~
n
.Ilt
.Il
n- t
lI.1lxn-t1I, where
a positive converging sequence with to = 0 and lim tn n~CO
that
Xn
converges to a limit x".
= r:
tn is
Show
152
5.
Applications to Numerical Analysis
IlllxnI s; ( Iltn)'Y (1IIlxn-lll)Y Ilt
5.6.
As in the previous exercise suppose that
5.7.
Show that the solution of
5.8.
Find the constant k in the previous problem and deduce that for Newton's method one has:
Zo -
n- I
~z~
IS given by Zn = 1 - Zn 1 - (1 - 2zo) 1/ 2 cotght cz"), with k appropriately chosen.
'I
X
*-
Xn
Zn+1
I s; -(1 2 f3y
=
2'
2zo)1/2 - -(J2 - " 1 - (J
where 1 - Zo - (1 - 2Z0) 1/ 2 (J = ---=--'-----=..:--:-= 1 - Zo (1- 2Z0 ) 1/ 2 '
+
5.9.
Obtain the result of Theorem 5.3.3 by considering the first integral with Zo = ~. (Hint: in this case the first order equation becomes Zn+1 = ~(1 + zJ.)
5.10. Consider the iterative method
= JL[AG(x n ) + (1- A)G(xk, Xk-I)] + (1- JL)Xk, A, JL E [0,1] and G, G are defined as in Theorem 5.3.5. Show Xk+1
where that JL and A can be chosen such that the convergence becomes faster. 5.11. Solve equation (5.5.1) directly and show the growth of the errors. 5.12. Obtain the relations (5.5.2), (5.5.3) and (5.5.4). 5.13. Show that the sweep method is equivalent to the Gaussian elimination of the problem Ay = g when the problem 5.5.1 is stated in vector form, see (5.4.16).
5.9. Notes The discussions of Sections 5.1 to 5.3 have been introduced only to show examples of application of difference equations to numerical analysis. More details can be found in the excellent books by Ortega and Rheinboldt [127] and Ostrowsky [132]. Theorem 5.3.3 is a simplified version of the original one, see [129], [127], while the estimates given before 5.3.4 (b) seem to be
5.9.
Notes
153
new. For other estimates of the error bounds for the Newton's method and Newton-like methods, see also [41], [109], [111], [179], [144]. For the effect of errors in the iterative methods, see Urabe [168]. A very large source of material on the problem discussed in Section 5.4 can be found in Gautschi's review paper [52]. Theorem 5.4.1 has been taken from Olver [123]. A detailed analysis, in a more general setting, of the Miller's algorithm is given in Zahar [180]. Applications to numerical methods for ODE's can be found in Cash's book [23], while a large exposition of applications as well as theoretical results are in Wimp's book [56], see also Mattheji [104]. The "sweep" method can be found in the Godunov and Rayabenki's book [59], where it is also presented for differential equations. An improvement of the Olver's algorithm can be found in [169]. An application of the sweep method to the problem of Section 5.4 can be found in Trigiante-Sivasundaram [167]. The contents of Section 5.6 are adapted from [89, 126, 127], while the material of Section 5.7 is taken from [91], see also [85, 87, 89].