Chapter 10 Numerical Solutions of Nonlinear Systems of Equations

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Chapter 10

Numerical Solutions of Nonlinear Systems of Equations

Hung-Yuan Fan (范洪源)

Department of Mathematics, National Taiwan Normal University, Taiwan

Spring 2016

Section 10.1

Fixed Points for Functions of Several

Variables

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Objective

To solve a system of nonlinear equations of the form











f₁(x1, x2,·, xn) = 0, f₂(x1, x2,·, xn) = 0, ...

f_n(x1, x2,·, xn) = 0,

(1)

where each f_i:Rⁿ→ R is a (nonlinear) function for i = 1, 2, . . . , n.

The unknown vector x = [x₁, x₂,· · · , xn]^T∈ Rⁿ is called a solution to the nonlinear system (1).

Vector-Valued Functions

Reformulation of the Nonlinear System

Consider a vector-valued function F :Rⁿ→ Rⁿ defined by F(x) = [f₁(x), f2(x),· · · , fn(x)]^T∈ Rⁿ ∀ x ∈ Rⁿ.

The system of nonlinear equations (1) can be represented as F(x) = 0, x = [x₁, x2,·, xn]^T∈ Rⁿ. (2) The functions f₁, f2,· · · , fnare called the coordinate

functions (坐標函數) of F.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Two Vector Norms (向量的範數)

Definition (常用的向量範數) Let v = [v₁, v₂,· · · , vn]^T∈ Rⁿ.

T he l₂

-norm (or Euclidean norm) of v is defined by

∥v∥2=√ v^Tv =

vu ut∑ⁿ

i=1

v²_i.

T he l_∞

-norm of v is defined by

1≤i≤n|vi|.

MATLAB Command:

∥v∥2, and norm(v, 'inf') is used for ∥v∥_∞.

Example 1, p. 631

Place the following nonlinear system







3x₁− cos(x2x₃)−¹₂ = 0

x²₁− 81(x2+ 0.1)²+ sin x3+ 1.06 = 0 e^−x1x2 + 20x3+^10π₃⁻³ = 0

(3)

in the form (2).

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Solution

Rewrite the nonlinear system (3) as

F(x₁, x₂, x₃)≡ [f1(x₁, x₂, x₃), f₂(x₁, x₂, x₃), f₃(x₁, x₂, x₃)]^T

= [0, 0, 0]^T= 0∈ R³, where the coordinate functions are defined by

f₁(x₁, x₂, x₃) = 3x₁− cos(x2x₃)−1 2,

f₂(x1, x2, x3) = x²₁− 81(x2+ 0.1)²+ sin x3+ 1.06, f₃(x₁, x₂, x₃) = e^−x1x2+ 20x₃+ 10π− 3

3 .

Fixed-Point Forms

As in Chap. 2, we shall transform the root-finding problem (2) into a fixed-point problem

x = G(x), x∈ D,

where G : D⊆ Rⁿ→ Rⁿ is some vector-valued function with domain

D ={[x1, x2,· · · , xn]^T| ai≤ xi≤ bi, i = 1, 2,· · · , n} (4) for some constants a₁, a₂,· · · , an and b₁, b₂,· · · , bn.

Fixed-Point Iteration (FPI) with an initial vector

provided that x^(k) ∈ D ∀ k ≥ 1.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Fixed Points in R

Def 10.5

The vector-valued function G : D⊆ Rⁿ→ Rⁿ has a fixed point at p∈ D ifG(p) = p.

Questions

Under what conditions does the sequence of vectors {x^(k)}^∞_k=0 generated by FPI converge to the unique fixed point p∈ D?

What is the error bound for the absolute error ∥x^(k)− p∥_∞? What is the rate of convergence for FPI?

We may write

G(x) = [g₁(x), g2(x),· · · , gn(x)]^T,

where each g_i is the ith component function of G for i = 1, 2, . . . , n.

Convergence Theorem of FPI

Thm 10.6

Let G be conti. on D withG(D)⊆ D, where the domain D is defined as in (4). Then

(1) G has at least one fixed point in D.

(2) If, in addition, ∃0 < K < 1s.t. each component function g_i has conti. partial derivatives with

∂gi(x)

∂xj

 ≤ K

n, whenever x∈ D,

for i, j = 1, 2, . . . , n, then with any x⁽⁰⁾∈ D conv. to the

unique fixed point p

∥x^(k)− p∥_∞≤ K^k

1− K∥x⁽¹⁾− x⁽⁰⁾∥_∞ ∀ k.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

How to check the continuity of G?

Thm 10.4 (分量函數的連續性)

Let g : D⊆ Rⁿ→ R be a function and x0 ∈ D. If ∃ δ > 0 and M > 0 s.t. the partial derivatives of g exist on N_δ(x0)∩ D with

∂g(x)

∂xj

 ≤ M ∀ x ∈ Nδ(x0)∩ D,

for j = 1, 2, . . . , n, theng is conti. at x₀. Continuity of G

G is conti. at x₀∈ D ⇐⇒ each component function gi is conti. at x₀ for i = 1, 2, . . . , n.

G is conti. on D⇐⇒ each gi is conti. on D for i = 1, 2, . . . , n.

Example 2, p. 633

(a) Place the nonlinear system in Example 1







3x₁− cos(x2x₃)−¹₂ = 0

x²₁− 81(x2+ 0.1)²+ sin x3+ 1.06 = 0 e^−x1x2 + 20x₃+^10π₃⁻³ = 0

in a fixed-point form x = G(x), x∈ D, and show that there is a unique sol. on

D ={[x1, x2, x3]^T| − 1 ≤ xi≤ 1, i = 1, 2, 3}.

(b) Perform the FPI withx⁽⁰⁾= [0.1, 0.1,−0.1]^T and the stopping criterion ∥x^(k)− x^(k−1)∥_∞< 10⁻⁵.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Solution of (a)

Solving the ith eq. of (3) for x_i (i = 1, 2, 3)⇒

x

3cos(x₂x₃) + 1

6 ≡g₁(x1, x2, x3)

x

9

√

x²₁+ sin x3+ 1.06− 0.1 ≡g₂(x1, x2, x3) (5)

x

20e^−x1x2 −10π− 3

60 ≡g₃(x₁, x₂, x₃).

So, define a vectored-valued function G : D→ R³ by

G(x₁, x2, x3) = [g₁(x1, x2, x3),g₂(x1, x2, x3),g₃(x1, x2, x3)]^T∈ R³ for any x = [x₁, x2, x3]^T∈ D. Now, consider the fixed-point form

x = G(x), x∈ D obtained from the original nonlinear system (3).

Solution of (a)–Conti’d

Fist, we shall claim thatG(D)⊆ D. It is easily seen from (5) that for any x∈ D, we have

|g1(x)| ≤ 1

3| cos(x2x₃)| +1

6 ≤ 0.50,

|g2(x)| =1 9

√

x²₁+ sin x3+ 1.06− 0.1

≤ 1 9

√(1)²+ sin(1) + 1.06− 0.1 < 0.09

|g3(x)| = 1

20e^−x1x2 +10π− 3 60

≤ 1

20e +10π− 3

60 < 0.61.

Hence, we know thatG(D)⊆ D.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Solution of (a)–Conti’d

Next, simple manipulation from Calculus gives that

∂g1

∂x2

= −x3

3 sin(x₂x₃), ∂g1

∂x3

= −x2

3 sin(x₂x₃), (6)

∂g2

∂x₁ = x₁

9√

x²₁+ sin x3+ 1.06, ∂g2

∂x₃ = cos x₃ 18√

x²₁+ sin x3+ 1.06, (7)

∂g₃

∂x1

= −x2

20 e^−x1x2, ∂g₃

∂x2

= −x1

20 e^−x1x2, ∂g₁

∂x1

= ∂g₂

∂x2

= ∂g₃

∂x3

= 0.

(8)

=⇒All first partial derivatives of g₁, g₂, g₃ are conti. on D!

Solution of (a)–Conti’d

Now, from (6) =⇒ ∂g1

∂x₂

 ≤ |x3|

3 · | sin(x2x₃)| ≤ sin 1

3 < 0.281, ∂g1

∂x₃

 < 0.281.

From (7), we see that ∂g₂

∂x1

 ≤ 1

9√

sin(−1) + 1.06 = 1 9√

0.218 < 0.238, ∂g₂

∂x₃

 ≤ 1

18√

sin(−1) + 1.06 = 1 18√

0.218 < 0.119, and furthermore, from (8), we also have

∂g₃

∂x1

 ≤ e

20 < 0.14, ∂g₃

∂x2

 ≤ e

20 < 0.14.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Solution of (a)–Conti’d

Thus, the partial derivatives of g₁, g₂, g₃ are bounded on D. It follows from Thm 10.4 thatG must be conti. on Dand

∂gi

∂xj

 ≤ 0.281 = K n = K

3 ∀ x ∈ D

for i, j = 1.2.3. So, the sufficient conditions of Thm 10.6 are satisfied with the constantK = (0.281)(3) = 0.843 < 1.

Conclusions

G has a unique fixed point p∈ D by Thm 10.6.

This fixed point p is one of the solutions to the original nonlinear system (3).

Solution of (b)–Numerical Results

Finally, perform the FPI

x^(k)= G(x^(k−1)), k = 1, 2, . . .

with x⁽⁰⁾ = [0.1, 0.1,−0.1]^T∈ D and∥x^(k)− x^(k−1)∥_∞< 10⁻⁵. Actual sol. p = [0.5, 0,^−π₆ ]^T ≈ [0.5, 0, −0.5235987757]^T.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

A Test for the Error Bound

With the computed sol. x⁽⁵⁾ and the actual fixed point p∈ D,

∥x⁽⁵⁾− p∥_∞≤ 2 × 10⁻⁸.

WithK = 0.843, the theoretical error bound would become

∥x⁽⁵⁾− p∥∞≤ (0.843)⁵

1−0.843(0.423) < 1.15.

The error bound in Thm 10.6 might be much larger than

the actual absolute error!

Accelerating Convergence (加速收斂性)

Basic Ideas

Use the latest estimates generated by the FPI x^(k)₁ , x^(k)₂ ,· · · , x^(k)_i−1

instead of x^(k−1)₁ , x^(k−1)₂ ,· · · , x^(k−1)_i−1 to compute the ith componentx^(k)_i .

This is the same as the idea of Gauss-Seidel method for solving linear systems. (See Chapter 7)

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Revisit Example 2

Reformulation as Gauss-Seidel Method

Consider the following Gauss-Seidel form for Example 2 x^(k)₁ = 1

3cos(x^(k₂ ⁻¹⁾x^(k₃⁻¹⁾) +1 6, x^(k)₂ = 1

9

√

(x^(k)₁ )²+ sin x^(k₃⁻¹⁾+ 1.06− 0.1, (9) x^(k)₃ = −1

20e^{−x(k)1 x(k)2} −10π− 3

60 , k = 1, 2, . . .

withx⁽⁰⁾= [0.1, 0.1,−0.1]^T∈ R³ and the same stopping criterion

∥x^(k)− x^(k−1)∥∞< 10⁻⁵.

Applying the iteration (9) with given initial vectorx⁽⁰⁾, the numerical results are shown in the following table. Note: In

general, this method does not always accelerate the

convergence! (不一定每次都能加速固定點迭代法!)

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Section 10.2

Newton’s Method

In One-Dimensional Case

Review of Newton’s Method

Newton’s method for solving a nonlinear equation of one variable

f(x) = 0, x∈ R

can be regarded as a fixed-point iteration with g(x) = x− 1

f^′(x) · f(x) ≡ x −ϕ(x)· f(x).

The quadratic convergence of Newton’s method is always expected if the initial guess is sufficiently close to a zero of f.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

In Multidimensional Case

Objectives

For solving a nonlinear system

F(x) = [f₁(x), f₂(x),· · · , fn(x)]^T= 0∈ Rⁿ, x∈ Rⁿ, try to develop a FPI with the vector-valued function

G(x) = x−A(x)⁻¹F(x)

≡ [g1(x), g2(x),· · · , gn(x)]^T∈ Rⁿ, x∈ Rⁿ, (10) assuming that A(x) = [a_ij(x)]∈ R^n×n is nonsingular at the

fixed point p of G.

Hopefully, the quadratic convergence can be achieved under reasonable conditions.

Thm 10.7 (FPI 二次收斂的充分條件) Let G(p) = p. Suppose that∃ δ > 0 with

(i) ^∂g_∂xⁱ

j is conti. on N_δ(p) for i, j = 1, 2, . . . , n;

(ii) _∂x^∂2gi

j∂xk is conti. on N_δ(p), and∃ M > 0 s.t.

∂²g_i(x)

∂xj∂x_k

 ≤ M ∀ x ∈ Nδ(p),

for i, j, k = 1, 2, . . . , n;

(iii) ^∂g_∂x^i(p)

j = 0 for i, j = 1, 2, . . . , n.

Then∃δˆ≤ δ s.t. the seq. {x^(k)}^∞_k=0 generated by FPI converges

quadratically to p for any x

∥x^(k)− p∥_∞≤ n²M

2 ∥x^(k−1)− p∥²_∞ ∀ k ≥ 1.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Derivation of the Matrix A(x)

WriteA(x)⁻¹= [bij(x)]∈ R^n×n. From (10)⇒

g_i(x) = xi−

∑n k=1

b_ik(x)f_k(x), i = 1, 2, . . . , n.

For each i, j = 1, 2, . . . , n, the first partial derivatives of g_i are

∂g_i(x)

∂xj

=









1− ∑ⁿ

k=1

(∂b_ik(x)

∂xj f_k(x) +b_ik(x)^∂f_∂x^k(x)_j )

, i = j,

−∑ⁿ

k=1

(∂bik(x)

∂xj f_k(x) +b_ik(x)^∂f_∂x^k(x)

j

)

, i̸= j.

(11)

Derivation of the Matrix A(x)–Conti’d

From condition (iii) of Thm 10.7 and (11), we immediately obtain

0 = ∂gi(p)

∂xj =









1− ∑ⁿ

k=1

b_ik(p)^∂f_∂x^k(p)_j , i = j,

−∑ⁿ

k=1

b_ik(p)^∂f_∂x^k(p)

j , i̸= j. (12) Define the Jacobian matrix J(x) = [^∂f_∂x^i(x)

j ]∈ R^n×n by

J(x) =









∂f1

∂x1(x) _∂x^∂f1

2(x) · · · _∂x^∂f1_n(x)

∂f2

∂x1(x) _∂x^∂f2

2(x) · · · _∂x^∂f2_n(x)

... ... ...

∂fn

∂x1(x) _∂x^∂fn

2(x) · · · _∂x^∂fn_n(x)







, x∈ Nδ(p).

It follows from (12) that A(p)⁻¹J(p) = I or A(p) = J(p).

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Newton’s Method

So, it is appropriate to choose A(x) = J(x)for x∈ Nδ(p).

Basic form of Newton’s method for nonlinear systems:

x^(k) = G(x^(k−1)) = x^(k−1)−A(x^(k−1))⁻¹F(x^(k−1))

= x^(k−1)−J(x^(k−1))⁻¹F(x^(k−1)), k = 1, 2, . . . , (13) where x⁽⁰⁾ ∈ Nˆδ(p) and J(x) is nonsingular on N_ˆδ(p) with 0 <δˆ≤ δ.

Quadratic convergence of Newton’s method is guaranteed

Some Comments on Newton’s Method (13)

We DO NOT computeJ(x^(k−1))⁻¹

explicitly in practical

In order to save the operation counts, we first solve the linear system

J(x^(k−1))y=−F(x^(k−1))

for the correction vectoryusing Gaussian Elimination

with Partial Pivoting, and then compute the next iterate via

x^(k)= x^(k−1)+y.

Floating-point operation counts ≈ O(²₃

n

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Pseudocode of Newton’s Method

To approx. the sol. of the nonlinear system F(x) = 0, x∈ Rⁿ. Algorithm 10.1: Newton’s Method for Systems

INPUT dim. n; initial x∈ Rⁿ; tol. TOL; max. no. of iter. N₀. OUTPUT an approx. sol. x to the nonlinear system.

Step 1 Set k = 1.

Step 2 While (k≤ N0) do Steps 3–7

Step 3 Compute F(x) and the Jacobian matrixJ(x).

Step 4 Solve the n× n linear systemJ(x)y=−F(x).

Step 5 Set x = x +y.

Step 6 If∥y∥ < TOL then OUTPUT(x); STOP.

Step 7 Set k = k + 1.

Step 8 OUTPUT(‘Maximum number of iterations exceeded’); STOP.

Example 1, p. 641 (See also Example 2 of Sec. 10.1) Apply Newton’s Method to solve the nonlinear system







3x₁− cos(x2x₃)−¹₂ = 0

x²₁− 81(x2+ 0.1)²+ sin x3+ 1.06 = 0 e^−x1x2+ 20x₃+ ^10π₃⁻³ = 0

withx⁽⁰⁾= [0.1, 0.1,−0.1]^T and∥x^(k)− x^(k−1)∥_∞< 10⁻⁵.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Numerical Results of Example 1

The Jacobian matrix J(x) is easily obtain from Calculus as

J(x₁, x2, x3) =



 3 x₃ ∼ (x2x₃) x₃sin(x₂x₃) 2x1 −162(x2+ 0.1) cos x₃

−x2e^−x1x2 −x1e^−x1x2 20



 .

Actual sol. p = [0.5, 0,^−π₆ ]^T≈ [0.5, 0, −0.5235987757]^T.

Section 10.3

Quasi-Newton Methods

(擬牛頓法)

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Newton’s Method v.s. Broyden’s Method (1/2)

For Each Iterate of Newton’s Method

At least n² scalar functional evaluations for the Jacobian matrix J(x^(k)) andn scalar functional evaluations for F(x^(k)).

Solving a linear system involving the Jacobian requires O(n³) operation counts.

Self-Correcting: it will generally correct for roundoff error

Quadratic convergenceoccurs if a good initial guess is given.

Newton’s Method v.s. Broyden’s Method (2/2)

For Each Iterate of Broyden’s Method

Only n scalar functional evaluations are required!

The amount of operation counts for solving the linear system is reduced to O(n²).

It is Not Self-Correcting with the successive iterations.

Only superlinear convergenceoccurs if a good initial guess is given, i.e., we have

klim→∞

∥x^(k+1)− p∥

∥x^(k)− p∥ = 0,

where p∈ Rⁿ is a solution of the nonlinear system F(x) = 0.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

About Broyden’s Method . . .

lt belongs to a class of least-change secant update

methods that produce algorithms called quasi-Newton.

The quasi-Newton methods replace the Jacobian matrix in Newton’s method with an approximate matrix that is

easily updated at each iteration.

References (參考文獻)

[Broy] C G. Broyden, A class of methods for solving nonlinear simultaneous equations, Math. Comp., 19(92), 577-593, 1965.

[DM] J. E. Dennis, Jr. and J. J. Moré, Quasi-Newton methods, motivation and theory, SIAM Rev., 19(1), 46–89, 1977.

Derivation of Broyden’s Method (1/2)

For an initial approx. x⁽⁰⁾∈ Rⁿ, compute the Jacobian matrix A₀ = J(x⁽⁰⁾)∈ R^n×n and the first iterate

x⁽¹⁾= x⁽⁰⁾−A₀⁻¹F(x⁽⁰⁾) as Newton’s method.

If we let

s₁= x⁽¹⁾− x⁽⁰⁾ and y₁= F(x⁽¹⁾)− F(x⁽⁰⁾), want to determine a matrixA₁ ≈ J(x⁽¹⁾)∈ R^n×n satisfying the quasi-Newton condition or secant condition

A₁(x⁽¹⁾− x⁽⁰⁾) = F(x⁽¹⁾)− F(x⁽⁰⁾) or A₁s₁ = y1. (14)

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Derivation of Broyden’s Method (2/2)

To determineA₁ uniquely, Broyden[Broy] imposed

A₁z = A₀z ∀ z ∈ Rⁿ with s^T₁z = 0 (15) on the secant codition (14).So, it follows from (14) and (15) that [DM]

A₁= A0+(y1− A0s₁)

∥s1∥²₂ · s^T1

and hence x⁽²⁾ = x⁽¹⁾−A₁⁻¹F(x⁽¹⁾).

In general, for k≥ 2, we have

A_k = Ak−1+(y_k− Ak−1s_k)

∥s_k∥²₂ ·s_k^T, (16) x^(k+1) = x^(k)−A_k⁻¹F(x^(k)),

wheres_k= x^(k)− x^(k−1)=−A⁻¹_k−1F(x^(k−1)) and y_k= F(x^(k))− F(x^(k−1)).

Remarks

From (16), we see that A_k is obtained from the previous A_k−1 plus an rank-1 updating matrix.

This technique is called the least-change secant updates.

In single-variable Newton’s method, may write f^′(x_k)≈ f(x_k)− f(xk−1)

x_k− xk−1 or f^′(x_k)(x_k−xk−1)≈ f(xk)−f(xk−1);

while we try to determine uniquelyA_k ≈ J(x^(k))s.t.

A_k(x^(k)− x^(k−1)) = F(x^(k))− F(x^(k−1)) in the multidimensional case.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

A Question

With the special structure of A_k, how to reduce the number of arithmetic calculations toO(n²) for solving the n× n linear system A_k⁻¹F(X^(k))?

Thm 10.8 (Sherman-Morrison Formula)

If A∈ R^n×n is nonsingular and x, y∈ Rⁿ are nonzero vectors with y^TA⁻¹x̸= −1, then A+xy^T is nonsingular and

(A+xy^T)⁻¹= A⁻¹−A⁻¹xy^TA⁻¹ 1 + y^TA⁻¹x.

Reformulation of A

For each k≥ 1, from (16) and Sherman-Morrison formula =⇒

A_k⁻¹ =

(A_k−1+(yk− Ak−1s_k)

∥sk∥²2

· s^Tk

)₋₁

= A⁻¹_k−1−A⁻¹_k−1(_y

k−Ak−1s_k

∥sk∥²2

)s^T_kA⁻¹_k−1

1 + s^T_kA⁻¹_k−1(

yk−Ak−1sk

∥sk∥²2

)

= A⁻¹_k−1− (A⁻¹_k−1y_k− sk)(s^T_kA⁻¹_k−1)

∥sk∥²₂+ s^T_kA⁻¹_k−1y_k− ∥sk∥²₂

= A⁻¹_k−1+(s_k− A⁻¹_k−1y_k)(s^T_kA⁻¹_k−1) s^T_kA⁻¹_k−1y_k

= A⁻¹_k−1+(sk−A⁻¹_k−1y_k)(s^T_kA⁻¹_k−1)

−s^T_k · (−A⁻¹_k−1y_k) . (17)

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Algorithm 10.2: Broyden’s Method

INPUT dim. n; initial x∈ Rⁿ; tol. TOL; max. no. of iter. N₀. OUTPUT an approx. sol. x of nonlinear system F(x) = 0.

Step 1 Set A₀ = J(x): the Jacobian matrix evaluated at x.

v = F(x). (Note: v = F(x⁽⁰⁾).) Step 2 Set A = A⁻¹₀ . (Use Gaussian elimination.)

Step 3 Set s =−Av; x = x + s; k = 1. (Note: s = s1, x = x⁽¹⁾.) Step 4 While (k≤ N0) do Steps 5–11.

Step 5 Set w = v; v = F(x); y = v− w. (Note: y = yk.) Step 6 Set z =−Ay. (Note: z = −A⁻¹k−1yk.)

Step 7 Set p =−s^Tz. (Note: p = s^T_kA⁻¹_k−1yk.)

Step 8 Set u^T= s^TA; A = A +¹_p(s + z)u^T. (Note: A = A⁻¹_k .) Step 9 Set s =−Av; x = x + s. (Note: s = −A⁻¹k F(x^(k)) and

x = x^(k+1).)

Step 10 If∥s∥ < TOL then OUTPUT(x); STOP.

Step 11 Set k = k + 1.

Step 12 OUTPUT(‘Maximum number of iterations exceeded’); STOP.

Example 1, p. 651 (See also Example 2 of Sec. 10.1) Use Broyden’s Method to solve the nonlinear system







3x₁− cos(x2x₃)−¹₂ = 0

x²₁− 81(x2+ 0.1)²+ sin x3+ 1.06 = 0 e^−x1x2+ 20x₃+ ^10π₃⁻³ = 0

withx⁽⁰⁾= [0.1, 0.1,−0.1]^T and∥x^(k)− x^(k−1)∥2 < 10⁻⁵.

. . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . .

. .

. . . . .

Numerical Results for Example 1

The superlinear convergence of Broyden’s method for Example 1 is demonstrated in the following table, and the computed solutions are less accurate than those computed by Newton’s method.

Actual sol. p = [0.5, 0,^−π₆ ]^T ≈ [0.5, 0, −0.5235987757]^T.

Thank you for your attention!