MATH 3795 Lecture 13. Numerical Solution of Nonlinear Equations in R .

MATH 3795
Lecture 13. Numerical Solution of Nonlinear
Equations in RN .
Dmitriy Leykekhman
Fall 2008
Goals
I
Learn about different methods for the solution of F (x) = 0, their
advantages and disadvantages.
I
Convergence rates.
I
MATLAB’s fsolve
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
1
Nonlinear Equations.
I
I
Goal: Given a function f : RN → RN we want to find x∗ ∈ RN
such that F (x∗ ) = 0.
Definition
A point x∗ with F (x∗ ) = 0 is called a root of F or zero of F .
I
We want to extend the methods from the last lecture, like Newton’s
method and Secant Method to find roots to vector valued problems.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
2
Convergence of Sequences.
Let {xk } be a sequence of vectors in RN and let k · k denote a vector
norm.
1. The sequence is called q-linearly convergent if there exists
c ∈ (0, 1) and k̂ ∈ N such that
kxk+1 − x∗ k ≤ ckxk − x∗ k
k ≥ k̂.
for all
2. The sequence is called q-superlinearly convergent if there exists a
sequence {ck } with ck > 0 and limk→∞ ck = 0 such that
kxk+1 − x∗ k ≤ ck kxk − x∗ k
or, equivalently, if
lim
k→∞
kxk+1 − x∗ k
= 0.
kxk − x∗ k
3. The sequence is called q-quadratically convergent to x∗ if
limk→∞ xk = x∗ and if there exists c > 0 and k̂ ∈ N such that
kxk+1 − x∗ k ≤ ckxk − x∗ k2
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
k ≥ k̂.
for all
Linear Least Squares
–
3
Taylor Expansion.
I
Given a differentiable function F , we define its Jacobian
 ∂F1

∂F1
∂x1 (x) . . .
∂xn (x)


..
..
F 0 (x) = JF (x) = 

.
.
∂Fn
∂Fn
∂x1 (x) . . .
∂xn (x)
I
If F is continuously differentiable around x0 , then
F (x) ≈ F (x0 ) + F 0 (x0 )(x − x0 )
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
4
Taylor Expansion.
Example
Let F : R2 → R2 be given by
F (x) =
x21 + x22 − 2
x1 −1
e
+ x32 − 2
Then its Jacobian is
0
2x1
ex1 −1
F (x) =
If x0 =
F (x) ≈
0
0
, then for x =
−2
e−1 − 2
+
0
e−1
x1
x2
0
0
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
2x2
3x22
close to
x1
x2
0
0
=
Linear Least Squares
–
−2
e−1 (1 + x1 ) − 2
5
Newton’s Method.
I
Suppose we are given an approximation x0 of a root x∗ of F .
I
Taylor approximation of F around x0 gives
F (x∗ ) = F (x0 + (x∗ − x0 )) ≈ F (x0 ) + F 0 (x0 )(x∗ − x0 ).
I
We use
M (x) = F (x0 ) + F 0 (x0 )(x − x0 )
I
as a model for F .
The root of M (x) i.e. the solution of linear system
F 0 (x0 )(x − x0 ) = −F (x0 )
is used as an approximation of the root of F .
I
Write the previous identity as
F 0 (x0 )s0 = −F (x0 ) step (correction)
x1 = x0 + s0 .
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
6
Newton’s Method.
Input: Initial values x(0), tolerance tol,
maximum number of iterations maxit
Output: approximation of the root
1 For k = 0:maxit do
2
Compute F’(x(k))s(k) = -F(x(k)).
3
Compute x(k+1) = x(k) + s(k).
4
Check for truncation
5 End
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
(LU-decomposition)
–
7
Newton’s Method.
Example
Consider
F (x) =
x21 + x22 − 2
x1 −1
e
+ x32 − 2
with starting point x0 = (1.5, 2)T and stopping criteria
kF (xk )k2 ≤ 1e − 10, the Newton’s method gives the computed solution
x = (1.0000, 1.0000)T and the history
k
kxk k2
kF (xk )k2
ksk k2
0 2.500000e + 000 8.750168e + 000 8.805454e − 001
1 1.665941e + 000 2.073196e + 000 3.234875e − 001
2 1.450739e + 000 4.127937e − 001 1.606253e − 001
3 1.423306e + 000 6.177196e − 002 2.206725e − 002
4 1.414386e + 000 1.401191e − 003 6.087256e − 004
5 1.414214e + 000 9.730653e − 007 3.964708e − 007
6 1.414214e + 000 4.415589e − 013 0.000000e + 000
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
8
Convergence of Newton’s Method.
Theorem
Let D ∈ RN be an open set and let F : D → RN be differentiable on D
with Lipschitz continuous derivative, i.e. let L > 0 be such thtat
kF 0 (y) − F 0 (x)k ≤ Lky − xk
∀x, y ∈ D.
If x∗ ∈ D is a root and if F 0 (x∗ ) is nonsingular, then there exists an
> 0 such that Newton’s method with starting point x0 with
kx0 − x∗ k < 0 generates iterates xk which converge to x∗ ,
lim xk = x∗ ,
k→∞
and which obey
kxk+1 − x∗ k ≤ ckxk − x∗ k2
for all k and some positive constant c.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
9
Convergence of Newton’s Method.
Theorem
Let D ∈ RN be an open set and let F : D → RN be differentiable on D
with Lipschitz continuous derivative, i.e. let L > 0 be such thtat
kF 0 (y) − F 0 (x)k ≤ Lky − xk
∀x, y ∈ D.
If x∗ ∈ D is a root and if F 0 (x∗ ) is nonsingular, then there exists an
> 0 such that Newton’s method with starting point x0 with
kx0 − x∗ k < 0 generates iterates xk which converge to x∗ ,
lim xk = x∗ ,
k→∞
and which obey
kxk+1 − x∗ k ≤ ckxk − x∗ k2
for all k and some positive constant c.
Newton’s method is locally q-quadratically convergent (under the
assumptions stated in the theorem).
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
9
Stopping Criteria.
We have discussed iterative methods which generate sequences xk with
(under suitable conditions) limk→∞ xk = x∗ . When do we stop the
iteration?
I
For a given tolerance tola > 0 we want to find xk such that
kxk − x∗ k < tola , stop if absolute error is small; or
kxk − x∗ k < tolr kx∗ k, stop if relative error is small.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
10
Stopping Criteria.
We have discussed iterative methods which generate sequences xk with
(under suitable conditions) limk→∞ xk = x∗ . When do we stop the
iteration?
I
For a given tolerance tola > 0 we want to find xk such that
kxk − x∗ k < tola , stop if absolute error is small; or
kxk − x∗ k < tolr kx∗ k, stop if relative error is small.
I
For some tk ∈ [0, 1]
kF (xk )k = kF (xk ) − F (x∗ )k
1
≥ k(F 0 (x∗ ))−1 k−1 kxk − x∗ k,
2
if xk is suff. close to x∗
Hence, if kF (xk )k < tolf , and if xk is sufficiently close to x∗ , then
kxk − x∗ k < 2tolf k(F 0 (x∗ ))−1 k.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
10
Stopping Criteria.
We have discussed iterative methods which generate sequences xk with
(under suitable conditions) limk→∞ xk = x∗ . When do we stop the
iteration?
I
For a given tolerance tola > 0 we want to find xk such that
kxk − x∗ k < tola , stop if absolute error is small; or
kxk − x∗ k < tolr kx∗ k, stop if relative error is small.
I
For some tk ∈ [0, 1]
kF (xk )k = kF (xk ) − F (x∗ )k
1
≥ k(F 0 (x∗ ))−1 k−1 kxk − x∗ k,
2
if xk is suff. close to x∗
Hence, if kF (xk )k < tolf , and if xk is sufficiently close to x∗ , then
kxk − x∗ k < 2tolf k(F 0 (x∗ ))−1 k.
I
Limit maximum number of iterations.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
10
Variations of Newton’s Method.
Recall Newton’s Method:
Input: Initial values x0 ,
tolerance tol, maximum number of iterations maxit
Output: approximation of the root
1. For k = 0, . . . , maxit do
2. Compute sk = −f (xk )/f 0 (xk ).
3. xk+1 = xk + sk .
4. Check for truncation
5. End
I
Requires the evaluation of derivatives and solution of a linear system.
I
Linear system solves are done using LU decomposition (cost 2/3n3
flops for each matrix factorization).
If Jacobian/derivative evaluations are expensive or difficult to
compute, or if the LU factorization of the derivative is expensive ,
the following variations are useful
I
I
Finite difference Newton method and Secant method.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
11
Finite Difference Newton Method.
Recall
Fi (x1 , . . . , xj−1 , xj + h, xj+1 , . . . , xn ) − Fi (x)
∂Fi (x)
= lim
.
h→0
∂xj
h
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
12
Finite Difference Newton Method.
Recall
Fi (x1 , . . . , xj−1 , xj + h, xj+1 , . . . , xn ) − Fi (x)
∂Fi (x)
= lim
.
h→0
∂xj
h
Thus, for all i = 1, . . . , n
 ∂F (x) 
1
∂xj



..
.
∂Fn (x)
∂xj

 = lim F (x1 , . . . , xj−1 , xj + h, xj+1 , . . . , xn ) − F (x) .
 h→0
h
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
12
Finite Difference Newton Method.
Recall
Fi (x1 , . . . , xj−1 , xj + h, xj+1 , . . . , xn ) − Fi (x)
∂Fi (x)
= lim
.
h→0
∂xj
h
Thus, for all i = 1, . . . , n
 ∂F (x) 
1
∂xj



..
.
∂Fn (x)
∂xj

 = lim F (x1 , . . . , xj−1 , xj + h, xj+1 , . . . , xn ) − F (x) .
 h→0
h
The idea is to replace
 ∂F (x) 
1
∂xj



..
.
∂Fn (x)
∂xj

 ≈ 1 F (x1 , . . . , xj−1 , xj + h, xj+1 , . . . , xn ) − F (x) .
 hj
Good choice for the step size hj = kxk.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
12
Finite Difference Newton Method.
Input: Initial values x0 ,
tolerance tol, maximum number of iterations maxit
Output: approximation of the root
1. For k = 0, . . . , maxit do
2. Compute finite difference approximation of F 0 (xk )
i.e., compute matrix B with columns B of F 0 (xk )
Bej =
1
F (x1 , . . . , xj−1 , xj + h, xj+1 , . . . , xn ) − F (x)
hj
3. Factor B
4. Solve Bsk = −F (xk )
5. Compute xk+1 = xk + sk
6. Check for truncation
7. End
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
13
Secant Method.
I
Recall that in the secant method we replace the derivative f 0 (xk+1 )
(xk+1 )
.
by f (xxkk)−f
−xk+1
I
If we set bk+1 =
f (xk )−f (xk+1 )
xk −xk+1
then bk+1 satisfies
bk+1 (xk+1 − xk ) = f (xk+1 ) − f (xk )
which is called the secant equation
I
The next iterate is given by
xk+2 = xk+1 + sk+1 ,
where sk+1 is the solution of bk+1 sk+1 = −f (xk+1 ).
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
14
Secant Method.
I
I
Try to extend this to the problem of finding a root of F : Rn → Rn .
Given two iterates xk , xk+1 ∈ Rn we try to find a nonsingular matrix
Bk+1 ∈ Rn×n which satisfies the so-called secant equation
Bk+1 (xk+1 − xk ) = F (xk+1 ) − F (xk ).
I
Then we compute the new iterate as follows
Solve
Bk+1 sk+1 = −F (xk+1 ),
xk+2 = xk+1 + sk+1 .
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
15
Secant Method.
There is a problem with this approach.
I
If n = 1 the the secant equation
bk+1 (xk+1 − xk ) = f (xk+1 ) − f (xk )
has a unique solution.
I
If n > 1, the we need to determine n2 entries of Bk+1 from n
equations
Bk+1 (xk+1 − xk ) = F (xk+1 ) − F (xk ).
There is no unique solution.
I
For example, for n = 2, xk+1 − xk = (1, 1)T and
F (xk+1 ) − F (xk ) = (1, 2), then the matrices
1 0
0 1
,
0 2
2 0
satisfy secant equation.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
16
Secant Method.
I
Therefore we chose Bk+1 ∈ Rn×n as the solution of
min
B(xk+1 −xk )=F (xk+1 )−F (xk )
I
kB − Bk kF .
Motivation: Bk+1 should satisfy the secant equation
B(xk+1 − xk ) = F (xk+1 ) − F (xk )
I
and Bk+1 should be as close to the old matrix Bk as possible to
preserve as much information contained in Bk as possible.
Common notation
sk = xk+1 − xk
and yk = F (xk+1 ) − F (xk ).
With these definition the previous minimization problem becomes
min
s:t: Bsk =yk
kB − Bk kF .
Unique solution
Bk+1 = Bk +
(yk − Bk sk )sTk
sTk sk
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Broyden update.
Linear Least Squares
–
17
Broyden’s Method.
Input: Initial values x0 ,
tolerance tol, maximum number of iterations maxit
Output: approximation of the root
1
2
3
4
5
6
For k = 0, . . . , maxit do
Solve Bk sk = −F (xk ) for sk
Set xk+1 = xk + sk .
Evaluate F (xk+1 )
Check for truncation
Set yk = F (xk+1 ) − F (xk )
7
8
Set Bk+1 = Bk +
End
(yk−Bk sk )sT
k
sT
k sk
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
.
Linear Least Squares
–
18
Broyden’s Method.
I
Note that due to the definition of yk and sk we have that
Bk+1 = Bk +
I
(yk − Bk sk )sTk
F (xk+1 )sTk
=
B
+
k
sTk sk
sTk sk
To start Broyden’s method we need an initial guess x0 for the root
x∗ and an initial matrix B0 ∈ Rn×n . In practice one often chooses
B0 = F 0 (x0 ) or
B0 = γI;
where γ is a suitable scalar. Other choices, for example finite
difference approximations to F 0 (x0 ) or choices based on the specific
structure of F 0 are also used.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
19
Broyden’s Method.
Example
Consider
F (x) =
x21 + x22 − 2
ex1 −1 + x32 − 2
with starting point
x0 =
1.5
2
and B0 = F 0 (x0 ),
and stopping criteria kF (xk )k2 ≤ 1e − 10.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
20
Broyden’s Method.
the Broyden’s method produces
k
0
1
2
3
4
5
6
7
8
9
10
kxk k2
2.500000e + 000
1.665941e + 000
1.476513e + 000
1.410326e + 000
1.417633e + 000
1.423860e + 000
1.415846e + 000
1.414375e + 000
1.414212e + 000
1.414214e + 000
1.414214e + 000
kF (xk )k2
8.750168e + 000
2.073196e + 000
8.734179e − 001
3.812507e − 001
1.586346e − 001
4.298504e − 002
4.681398e − 003
6.074087e − 004
4.051447e − 006
2.724111e − 008
1.182169e − 011
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
ksk k2
8.805454e − 001
1.922038e − 001
1.321894e − 001
1.555213e − 001
9.620188e − 002
1.043037e − 002
2.583147e − 003
6.288185e − 004
1.805771e − 006
1.154246e − 008
0.000000e + 000
Linear Least Squares
–
21
MATLAB’s fsolve.
Similar to a build-in function fzero, Optimization Toolbox has a function
fsolve that tries to solve a system of nonlinear equations
F (x) = 0.
Warning: You need this Toolbox in order to use this function.
Syntax:
x = fsolve(fun,x0)
x = fsolve(fun,x0,options)
x = fsolve(problem)
[x,fval] = fsolve(fun,x0)
[x,fval,exitflag] = fsolve(...)
[x,fval,exitflag,output] = fsolve(...)
[x,fval,exitflag,output,jacobian] = fsolve(...)
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
22
MATLAB’s fsolve.
x = fsolve(fun,x0)
starts at x0 and tries to solve the equations described in fun.
[x,fval] = fsolve(fun,x0)
returns the value of the objective function fun at the solution x.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
23
MATLAB’s fsolve.
Thus in our example with
x21 + x22 − 2
F (x) =
ex1 −1 + x32 − 2
and x0 =
[x,fval] = fsolve(’ex1’,[1.5 2])
produces
x =
fval =
0.999999999999919
1.000000000000164
1.0e-012 *
0.165201186064223
0.410338429901458
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
24
1.5
2
MATLAB’s fsolve.
Example
Let’s say you want to solve
2x1 − x2 = e−x1
2
−x1 + 2x2 = e−x .
In other words we want to solve nonlinear system
2x1 − x2 − e−x1 = 0
2
−x1 + 2x2 − e−x = 0.
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
25
MATLAB’s fsolve.
First we create a function in a m-file.
function varargout =
%
%
% [F]
= ex3(x)
% [F,Jac] = ex3(x)
%
%
ex3(x)
returns the function value in F
returns the function value in F and
the Jacobian in Jac
% return the function value
varargout{1} = [ 2*x(1)-x(2)-exp(-x(1));
-x(1) + 2*x(2) - exp(-x(2))];
if( nargout > 1 )
% return the Jacobian as the second argument
varargout{2} = [ 2+exp(-x(1))
-1;
-1
2+exp(-x(2))];
end
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
26
MATLAB’s fsolve.
Now calling
[x,fval] = fsolve(’ex3’,[-5 -5])
produces
x =
0.567143031397357
0.567143031397357
fval =
1.0e-006 *
-0.405909605705190
-0.405909605705190
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
27
MATLAB’s fsolve.
We got 6 digits of accuracy. If we need more we can use
options = optimset(’TolFun’,1e-10)
then calling
[x,fval] = fsolve(’ex3’,[-5 -5],options)
produces
x =
0.567143290409772
0.567143290409772
fval =
1.0e-013 *
-0.179856129989275
-0.179856129989275
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
28
MATLAB’s fsolve.
You could learn more about fsolve by typing
help fsolve
or looking at the function on mathworks website fsolve
Similarly you learn there more about the options by typing
help optimset
or just
optimset
D. Leykekhman - MATH 3795 Introduction to Computational Mathematics
Linear Least Squares
–
29