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