<< . .

. 41
( : 95)



. . >>

for the Dekker-Brent algorithm (for the meaning of ξ, see Example 6.5),
while for all the other schemes the search interval has been taken equal to
[0, 0.5]. The stopping tolerance toll has been set to 10’10 . The obtained
results are reported in Table 6.10 where nit and f (nit) denote respectively
the number of iterations needed by the method to converge and the value
of f at the computed solution.
Notice the extremely slow convergence of the Regula Falsi method, due
to the fact that the value v (k ) always coincides with the right end-point
v = 0.5 and the function f around v — has derivative very close to zero. An
analogous interpretation holds for the chord method.

f (nit) f (nit)
Method Method
nit nit
’1.12 · 10’15 1.09 · 10’14
bisection 33 Dekker-Brent 11
’9.77 · 10’11 2.7 · 10’20
Regula Falsi 225 secant 11
’9.80 · 10’14 ’1.35 · 10’20
chord 186 Newton™s 8
TABLE 6.10. Convergence of the methods for the approximation of the root of
equation (6.42)




6.8 Exercises
1. Derive geometrically the sequence of the ¬rst iterates computed by bisec-
tion, Regula Falsi, secant and Newton™s methods in the approximation of
the zero of the function f (x) = x2 ’ 2 in the interval [1, 3].
2. Let f be a continuous function that is m-times di¬erentiable (m ≥ 1), such
that f (±) = . . . = f (m’1) (±) = 0 and f (m) (±) = 0. Prove (6.22) and check
that the modi¬ed Newton method (6.23) has order of convergence equal to
2.
[Hint: let f (x) = (x ’ ±)m h(x), h being a function such that h(±) = 0].
3. Let f (x) = cos2 (2x) ’ x2 be the function in the interval 0 ¤ x ¤ 1.5
examined in Example 6.4. Having ¬xed a tolerance µ = 10’10 on the abso-
lute error, determine experimentally the subintervals for which Newton™s
method is convergent to the zero ± 0.5149.
[Solution: for 0 < x(0) ¤ 0.02, 0.94 ¤ x(0) ¤ 1.13 and 1.476 ¤ x(0) ¤ 1.5,
the method converges to the solution ’±. For any other value of x(0) in
[0, 1.5], the method converges to ±].
280 6. Root¬nding for Nonlinear Equations

4. Check the following properties:
(a) 0 < φ (±) < 1: monotone convergence, that is, the error x(k) ’ ±
maintains a constant sign as k varies;
(b) ’1 < φ (±) < 0: oscillatory convergence that is, x(k) ’ ± changes sign
as k varies;
(c) |φ (±)| > 1: divergence. More precisely, if φ (±) > 1, the sequence
is monotonically diverging, while for φ (±) < ’1 it diverges with
oscillatory sign.
5. Consider for k ≥ 0 the ¬xed-point method, known as Ste¬ensen™s method
f (x(k) + f (x(k) )) ’ f (x(k) )
f (x(k) )
x(k+1) = x(k) ’ •(x(k) ) =
, ,
•(x(k) ) f (x(k) )
and prove that it is a second-order method. Implement the Ste¬ensen
method in a MATLAB code and employ it to approximate the root of
the nonlinear equation e’x ’ sin(x) = 0.
6. Analyze the convergence of the ¬xed-point method x(k+1) = φj (x(k) ) for
computing the zeros ±1 = ’1 and ±2 = 2 of the function f (x) = x2 ’ x ’ 2,
when the following iteration functions are used: φ1 (x) = x2 ’ 2, φ2 (x) =
√ √
2 + x φ3 (x) = ’ 2 + x and φ4 (x) = 1 + 2/x, x = 0.
[Solution: the method is non convergent with φ1 , it converges only to ±2 ,
with φ2 and φ4 , while it converges only to ±1 with φ3 ].
7. For the approximation of the zeros of the function f (x) = (2x2 ’ 3x ’
2)/(x ’ 1), consider the following ¬xed-point methods:
(1) x(k+1) = g(x(k) ), where g(x) = (3x2 ’ 4x ’ 2)/(x ’ 1);
(2) x(k+1) = h(x(k) ), where h(x) = x ’ 2 + x/(x ’ 1).
Analyze the convergence properties of the two methods and determine
in particular their order. Check the behavior of the two schemes using
Program 51 and provide, for the second method, an experimental estimate
of the interval such that if x(0) is chosen in the interval then the method
converges to ± = 2.
[Solution: zeros: ±1 = ’1/2 and ±2 = 2. Method (1) is not convergent,
while (2) can approximate only ±2 and is second-order. Convergence holds
for any x(0) > 1].
8. Propose at least two ¬xed-point methods for approximating the root ±
0.5885 of equation e’x ’ sin(x) = 0 and analyze their convergence.
9. Using Descartes™s rule of signs, determine the number of real roots of the
polynomials p6 (x) = x6 ’ x ’ 1 and p4 (x) = x4 ’ x3 ’ x2 + x ’ 1.
[Solution: both p6 and p4 have one negative and one positive real root].

10. Let g : R ’ R be de¬ned as g(x) = 1 + x2 . Show that the iterates of
Newton™s method for the equation g (x) = 0 satisfy the following proper-
ties:
|x(0) | < 1 ’ g(x(k+1) ) < g(x(k) ), k ≥ 0, lim x(k) = 0,
(a)
k’∞
|x(0) | > 1 ’ g(x(k+1) ) > g(x(k) ), k ≥ 0, lim |x(k) | = +∞.
(b)
k’∞
7
Nonlinear Systems and Numerical
Optimization




In this chapter we address the numerical solution of systems of nonlinear
equations and the minimization of a function of several variables.
The ¬rst problem generalizes to the n-dimensional case the search for
the zeros of a function, which was considered in Chapter 6, and can be
formulated as follows: given F : Rn ’ Rn ,

¬nd x— ∈ Rn such that F(x— ) = 0. (7.1)

Problem (7.1) will be solved by extending to several dimensions some of
the schemes that have been proposed in Chapter 6.
The basic formulation of the second problem reads: given f : Rn ’ R,
called an objective function,

minimize f (x) in Rn , (7.2)

and is called an unconstrained optimization problem.
A typical example consists of determining the optimal allocation of n
resources, x1 , x2 , . . . , xn , in competition with each other and ruled by a
speci¬c law. Generally, such resources are not unlimited; this circumstance,
from a mathematical standpoint, amounts to requiring that the minimizer
of the objective function lies within a subset „¦ ‚ Rn , and, possibly, that
some equality or inequality constraints must be satis¬ed.
When these constraints exist the optimization problem is called con-
strained and can be formulated as follows: given the objective function f ,

minimize f (x) in „¦ ‚ Rn . (7.3)
282 7. Nonlinear Systems and Numerical Optimization

Remarkable instances of (7.3) are those in which „¦ is characterized by con-
ditions like h(x) = 0 (equality constraints) or h(x) ¤ 0 (inequality con-
straints), where h : Rn ’ Rm , with m ¤ n, is a given function, called cost
function, and the condition h(x) ¤ 0 means hi (x) ¤ 0, for i = 1, . . . , m.
If the function h is continuous and „¦ is connected, problem (7.3) is
usually referred to as a nonlinear programming problem. Notable examples
in this area are:

convex programming if f is a convex function and h has convex compo-
nents (see (7.21));

linear programming if f and h are linear;

quadratic programming if f is quadratic and h is linear.


Problems (7.1) and (7.2) are strictly related to one another. Indeed, if we
denote by Fi the components of F, then a point x— , a solution of (7.1),
n
is a minimizer of the function f (x) = i=1 Fi2 (x). Conversely, assuming
that f is di¬erentiable and setting the partial derivatives of f equal to
zero at a point x— at which f is minimum leads to a system of nonlinear
equations. Thus, any system of nonlinear equations can be associated with
a suitable minimization problem, and vice versa. We shall take advantage
of this observation when devising e¬cient numerical methods.



7.1 Solution of Systems of Nonlinear Equations
Before considering problem (7.1), let us set some notation which will be
used throughout the chapter.
For k ≥ 0, we denote by C k (D) the set of k-continuously di¬erentiable
functions from D to Rn , where D ⊆ Rn is a set that will be made precise
from time to time. We shall always assume that F ∈ C 1 (D), i.e., F : Rn ’
Rn is a continuously di¬erentiable function on D.
We denote also by JF (x) the Jacobian matrix associated with F and
evaluated at the point x = (x1 , . . . , xn )T of Rn , de¬ned as

‚Fi
(JF (x))ij = (x), i, j = 1, . . . , n.
‚xj

Given any vector norm · , we shall henceforth denote the sphere of radius
R with center x— by

B(x— ; R) = {y ∈ Rn : y ’ x— < R} .
7.1 Solution of Systems of Nonlinear Equations 283

7.1.1 Newton™s Method and Its Variants
An immediate extension to the vector case of Newton™s method (6.16) for
scalar equations can be formulated as follows:

given x(0) ∈ Rn , for k = 0, 1, . . . , until convergence:

JF (x(k) )δx(k) = ’F(x(k) );
solve
(7.4)
x(k+1) = x(k) + δx(k) .
set

Thus, at each step k the solution of a linear system with matrix JF (x(k) )
is required.

Example 7.1 Consider the nonlinear system
±
 ex2 +x2 ’ 1 = 0,
1 2


 ex2 ’x2 ’ 1
1 2 = 0,
2 2
which admits the unique solution x— = 0. In this case, F(x) = (ex1 +x2 ’
2 2
1, ex1 ’x2 ’ 1). Running Program 57, leads to convergence in 15 iterations to the
pair (0.61 · 10’5 , 0.61 · 10’5 )T , starting from the initial datum x(0) = (0.1, 0.1)T ,
thus demonstrating a fairly rapid convergence rate. The results, however, dra-
matically change as the choice of the initial guess is varied. For instance, picking
up x(0) = (10, 10)T , 220 iterations are needed to obtain a solution comparable to
the previous one, while, starting from x(0) = (20, 20)T , Newton™s method fails to

converge.

The previous example points out the high sensitivity of Newton™s method
on the choice of the initial datum x(0) , as con¬rmed by the following local
convergence result.

Theorem 7.1 Let F : Rn ’ Rn be a C 1 function in a convex open set
D of Rn that contains x— . Suppose that J’1 (x— ) exists and that there exist
F
positive constants R, C and L, such that J’1 (x— ) ¤ C and
F

∀x, y ∈ B(x— ; R),
JF (x) ’ JF (y) ¤ L x ’ y

having denoted by the same symbol · two consistent vector and matrix
norms. Then, there exists r > 0 such that, for any x(0) ∈ B(x— ; r), the
sequence (7.4) is uniquely de¬ned and converges to x— with

x(k+1) ’ x— ¤ CL x(k) ’ x— 2 . (7.5)
Proof. Proceeding by induction on k, let us check (7.5) and, moreover, that
x(k+1) ∈ B(x— ; r), where r = min(R, 1/(2CL)). First, we prove that for any
x(0) ∈ B(x— ; r), the inverse matrix J’1 (x(0) ) exists. Indeed
F

1
J’1 (x— )[JF (x(0) ) ’ JF (x— )] ¤ J’1 (x— ) JF (x(0) ) ’ JF (x— ) ¤ CLr ¤ ,
F F
2
284 7. Nonlinear Systems and Numerical Optimization

and thus, thanks to Theorem 1.5, we can conclude that J’1 (x(0) ) exists, since
F


J’1 (x— )
J’1 (x(0) ) ¤ ¤ 2 J’1 (x— ) ¤ 2C.
F
F F
1 ’ J’1 (x— )[JF (x(0) ) ’ JF (x— )]
F


As a consequence, x(1) is well de¬ned and

x(1) ’ x— = x(0) ’ x— ’ J’1 (x(0) )[F(x(0) ) ’ F(x— )].
F


Factoring out J’1 (x(0) ) on the right hand side and passing to the norms, we get
F


¤ J’1 (x(0) )
x(1) ’ x— F(x— ) ’ F(x(0) ) ’ JF (x(0) )[x— ’ x(0) ]
F

L—
¤ 2C x ’ x(0) 2
2

where the remainder of Taylor™s series of F has been used. The previous relation
proves (7.5) in the case k = 0; moreover, since x(0) ∈ B(x— ; r), we have x— ’
1—
x(0) ¤ 1/(2CL), from which x(1) ’ x— ¤ x ’ x(0) .
2
This ensures that x(1) ∈ B(x— ; r).
By a similar proof, one can check that, should (7.5) be true for a certain k,
then the same inequality would follow also for k + 1 in place of k. This proves
3
the theorem.

Theorem 7.1 thus con¬rms that Newton™s method is quadratically conver-
gent only if x(0) is su¬ciently close to the solution x— and if the Jacobian
matrix is nonsingular. Moreover, it is worth noting that the computational
e¬ort needed to solve the linear system (7.4) can be excessively high as n
gets large. Also, JF (x(k) ) could be ill-conditioned, which makes it quite di¬-
cult to obtain an accurate solution. For these reasons, several modi¬cations
to Newton™s method have been proposed, which will be brie¬‚y considered
in the later sections, referring to the specialized literature for further details
(see [OR70], [DS83], [Erh97], [BS90] and the references therein).


7.1.2 Modi¬ed Newton™s Methods
Several modi¬cations of Newton™s method have been proposed in order
to reduce its cost when the computed solution is su¬ciently close to x— .
Further variants, that are globally convergent, will be introduced for the
solution of the minimization problem (7.2).

1. Cyclic updating of the Jacobian matrix
An e¬cient alternative to method (7.4) consists of keeping the Jacobian
matrix (more precisely, its factorization) unchanged for a certain number,
say p ≥ 2, of steps. Generally, a deterioration of convergence rate is accom-
panied by a gain in computational e¬ciency.
7.1 Solution of Systems of Nonlinear Equations 285

Program 57 implements Newton™s method in the case in which the LU
factorization of the Jacobian matrix is updated once every p steps. The pro-
grams used to solve the triangular systems have been described in Chapter
3.
Here and in later codings in this chapter, we denote by x0 the initial
vector, by F and J the variables containing the functional expressions of F
and of its Jacobian matrix JF , respectively. The parameters toll and nmax
represent the stopping tolerance in the convergence of the iterative process
and the maximum admissible number of iterations, respectively. In output,
the vector x contains the approximation to the searched zero of F, while
nit denotes the number of iterations necessary to converge.


Program 57 - newtonxsys : Newton™s method for nonlinear systems

function [x, nit] = newtonsys(F, J, x0, toll, nmax, p)
[n,m]=size(F); nit=0; Fxn=zeros(n,1); x=x0; err=toll+1;
for i=1:n, for j=1:n, Jxn(i,j)=eval(J((i-1)*n+j,:)); end; end
[L,U,P]=lu(Jxn); step=0;
while err > toll
if step == p
step = 0;
for i=1:n;
Fxn(i)=eval(F(i,:));
for j=1:n; Jxn(i,j)=eval(J((i-1)*n+j,:)); end
end
[L,U,P]=lu(Jxn);
else
for i=1:n, Fxn(i)=eval(F(i,:)); end
end
nit=nit+1; step=step+1; Fxn=-P*Fxn; y=forward col(L,Fxn);
deltax=backward col(U,y); x = x + deltax; err=norm(deltax);
if nit > nmax
disp(™ Fails to converge within maximum number of iterations ™);
break
end
end



2. Inexact solution of the linear systems
Another possibility consists of solving the linear system (7.4) by an iter-
ative method where the maximum number of admissible iterations is ¬xed
a priori. The resulting schemes are identi¬ed as Newton-Jacobi, Newton-
SOR or Newton-Krylov methods, according to the iterative process that is
used for the linear system (see [BS90], [Kel99]). Here, we limit ourselves to
describing the Newton-SOR method.
286 7. Nonlinear Systems and Numerical Optimization

In analogy with what was done in Section 4.2.1, let us decompose the
Jacobian matrix at step k as

JF (x(k) ) = Dk ’ Ek ’ Fk (7.6)

where Dk = D(x(k) ), ’Ek = ’E(x(k) ) and ’Fk = ’F(x(k) ), the diagonal
part and the lower and upper triangular portions of the matrix JF (x(k) ),
respectively. We suppose also that Dk is nonsingular. The SOR method for
(k)
solving the linear system in (7.4) is organized as follows: setting δx0 = 0,
solve

δx(k) = Mk δxr’1 ’ ωk (Dk ’ ωk Ek )’1 F(x(k) ),
(k)
r = 1, 2, . . . , (7.7)
r

where Mk is the iteration matrix of SOR method
’1
Mk = [Dk ’ ωk Ek ] [(1 ’ ωk )Dk + ωk Fk ] ,

and ωk is a positive relaxation parameter whose optimal value can rarely
be determined a priori. Assume that only r = m steps of the method are
(k) (k)
carried out. Recalling that δxr = xr ’ x(k) and still denoting by x(k+1)
the approximate solution computed after m steps, we ¬nd that this latter
can be written as (see Exercise 1)
’1
x(k+1) = x(k) ’ ωk Mm’1 + · · · + I (Dk ’ ωk Ek ) F(x(k) ). (7.8)
k

This method is thus a composite iteration, in which at each step k, starting
from x(k) , m steps of the SOR method are carried out to solve approxi-
mately system (7.4).
The integer m, as well as ωk , can depend on the iteration index k; the
simplest choice amounts to performing, at each Newton™s step, only one
iteration of the SOR method, thus obtaining for r = 1 from (7.7) the one-
step Newton-SOR method
’1
x(k+1) = x(k) ’ ωk (Dk ’ ωk Ek ) F(x(k) ).

In a similar way, the preconditioned Newton-Richardson method with ma-
trix Pk , if truncated at the m-th iteration, is

P’1 F(x(k) ),
x(k+1) = x(k) ’ I + Mk + . . . + Mk
m’1
k

where Pk is the preconditioner of JF and

Mk = P’1 Nk , Nk = Pk ’ JF (x(k) ).
k

<< . .

. 41
( : 95)



. . >>