Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Lecture on Numerical Analysis
Dr.-Ing. Michael Dumbser
M. Dumbser
1 / 23
08 / 11 / 2007
Analisi Numerica
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Numerical Solution of Nonlinear Equations
Task: compute approximately the solution x of the nonlinear algebraic equation
f ( x)  0
(1)
If the problem is not in form (1), then it can always be cast into form (1) by putting all
terms on the left hand side.
Examples:
exp( x )  sin( x )  x 2
exp( x )  sin( x )  x 2  0
f ( x )  exp( x )  sin( x )  x 2
g ( x )  h ( x )  e( x )
f ( x )  g ( x )  h ( x )  e( x )
M. Dumbser
2 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Numerical Solution of Nonlinear Equations
Task: compute approximately the solution x of the nonlinear algebraic equation
f ( x)  0
(1)
Two problems present even at the analytical level:
• For general functions f(x), there may be no solution at all.
Example:
f ( x )  sin( x )  2  0
• For general functions f(x), there may be more than one solution.
Example:
f ( x)  x 2  1  0
M. Dumbser
3 / 23
x1  1,
x2  1
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Regula Falsi Method
f
x
x
0
2
x4
x
x
3
1
x
(I) Before we can start using this method, we must guess two points x0 and x1 where the function f(x)
has opposite sign, i.e.
f ( x0 )  f ( x1 )  0
(II) Now connect the two points (x0, f(x0)) and (x1, f(x1)) with a straight line and compute the intersection
of this line with the x-axis, denoted by xi. The point xi is given by
xi 
f 0 x1  x0 f1
f 0  f1
(III) Compute the sign of f(x2) and compare against the signs of f(x0) and f(x1). Then update the interval
boundaries according to the result, i.e. set
x0 : xi
M. Dumbser
4 / 23
if
f ( x0 )  f ( xi )  0
x1 : xi
if
f ( x1 )  f ( xi )  0
If |f(xi)| < tolerance then STOP (in Maple: break), otherwise go back to step (II).
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Regula Falsi Method
f
x
x
0
2
x4
x
x
3
1
x
Advantages: - it is guaranteed that the method finds a solution of f(x)=0 because we already know that it
the solution is contained in the initial interval [x0;x1].
- The method is very robust and reliable
Disadvantage: - it can not be applied if we do not have an initial guess of the interval [x0;x1].
M. Dumbser
5 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Secant Method
f
x
x x
x
x
2
0
1
3
• Variant of the regula falsi method: we start from any two points x0 and x1
• Advantage: no initial interval must be guessed
• Disadvantage: may not always converge towards a solution !
M. Dumbser
6 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
The Secant Method
f
x
x
x
0
1
• Variant of the regula falsi method: we start from any two points x0 and x1
• Advantage: no initial interval must be guessed
• Disadvantage: may not always converge towards a solution !
M. Dumbser
7 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Newton‘s Method
f
x
x
1
x2
x
0
• Differential form of the secant method. Instead of using the secant, Newton proposes to use the
tangent!
• Advantages: - no initial interval must be guessed
- very fast
• Disadvantage: may not always converge towards a solution !
M. Dumbser
8 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Newton‘s Method
f
x
x
1
x2
x
0
• Differential form of the secant method. Instead of using the secant, Newton proposes to use the
tangent!
f
f ( xi )   x  xi   0
x
f
x   f ( xi ) /
x
xi 1  xi  x
M. Dumbser
9 / 23
(1)
(2)
x : x  xi
If |f(xi)| < tolerance then STOP (in Maple/MATLAB: break) or
If | x | < tolerance then STOP
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 1
Write a Maple worksheet and a MATLAB script that solve the nonlinear equation
f ( x )  exp( x)  sin( x )  x 2  0
using the regula falsi method. Choose an appropriate starting interval.
Exercise 2
Now write a MATLAB script that solves the problem of exercise 1 using the Newton
method. The user should have the possibility to choose the starting guess, the
maximum number of iterations as well as the tolerance.
M. Dumbser
10 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Globally Convergent Newton Methods
From the previous Example 4 we find that even for strictly monotone functions f(x) Newton‘s
method may not always converge from arbitrary starting points x0, but only from starting
points that are sufficiently close to the solution of the problem. However, if the starting point
is sufficiently close, the method converges. We therefore say the original Newton method is
only locally convergent.
The local convergence behaviour of Newton‘s method and the possible failure to converge
for initial guesses far away from the solution can be mainly explained by too large step sizes
x:
M. Dumbser
11 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Globally Convergent Newton Methods
As seen in the sketch of the previous slide, the absolute value of the function f(x) slightly
increased during one Newton step. To avoid this problem, we therefore add the following
requirement: The function |f(x)| must always decrease during one Newton step, i.e.:
f ( xi  xi )  f ( xi )
If this is not fulfilled, the step size is sucessively reduced until the criterion is satisfied. The
globally convergent Newton method then reads as follows:
DO i=0,1,…
f
x   f ( xi ) /
x
 1
WHILE f ( xi  xi )  f ( xi )
   
END
xi 1  xi    x
M. Dumbser
12 / 23
UNTIL
f ( xi )  tol
0  1
Free parameter.
Choose e.g.  = 0.5
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Globally Convergent Newton Methods
Exercise 5
Write a modified MATLAB script GlobalNewton.m to solve the nonlinear equation
f ( x )  arctan( x  1)  0
with the globally convergent Newton method.
1)
Use x0 = 2 as initial guess.
2)
Use x0 = 2.5 as initial guess. What do you realize?
3)
Use x0 = 106 as initial guess. Does the method converge?
M. Dumbser
13 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Computation of the Roots of a Polynomial
We now apply Newton‘s method to a special case of nonlinear algebraic equation where the
function f(x) is a polynomial of degree n:
n
f ( x )  pn ( x )   ai x i  0
i 0
The values x which satisfy eqn. (1) are called the roots of the polynomial pn(x).
We suppose in the following that all the roots of the polynomial are distinct and real.
Suppose x = x is a root of pn(x), i.e.
pn (x )  0
then we have
pn ( x )  ( x  x ) pn 1 ( x )
If we knew all the roots xi we would have the following expression of the polynomial:
pn ( x)  ( x  x1 )( x  x2 )...( x  xn )
M. Dumbser
14 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Computation of the Roots of a Polynomial
A simple and straightforward application of Newton‘s method now consists in finding a first
root x1 (typically the largest one) of the polynomial pn and then compute the largest root x2 of
the polynomial pn-1 obtained by dividing the polynomial pn by ( x – x ):
pn 1 ( x ) 
pn ( x )
x  x1
After having found the root x2 of pn-1 we divide again
pn 2 ( x ) 
pn 1 ( x )
pn ( x )

x  x2
( x  x1 )( x  x 2 )
and compute the next root x3 and so on, until all the n roots of pn have been found.
In general, we have:
pn ( x )
pn  j ( x ) 
( x  x1 )...( x  x j )
j


pn ( x )
1
x pn ( x )
pn  j ( x ) 


x
( x  x1 )...( x  x j ) ( x  x1 )...( x  x j ) i 1 x  xi
M. Dumbser
15 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Computation of the Roots of a Polynomial
To express the polynomial pn-j(x) and its derivative, two simple MATLAB loops defining the
auxiliary variables Dj and Sj are necessary:
pn ( x )
Dj


p ( x)
x pn ( x )
pn j ( x) 
 n
Sj
x
Dj
Dj
pn  j ( x ) 
with
Dj  1
Sj  0
FOR k  [1 : j]
D j  D j * ( x  xk )
FOR k  [1 : j]
S j  S j  1 /( x  xk )
END
END
We then obtain the so-called Newton-Maehly method:
x   f ( x ) /
M. Dumbser
16 / 23
f

  pn j ( xi ) / pn j 
x
x

x
 pn
pn  pn S j
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 3
(1) Write a MATLAB script PolyRoot.m that finds all the roots of a general polynomial
of degree n.
The user should be able to specify the coefficients a of the polynomial. Then, two
MATLAB function files, called Polynomial.m and DPolynomial.m compute the
polynomial and its first derivative in point x using the coefficients a.
Further user inputs required are tol, NMAX and x0.
(2) Apply your general MATLAB script to the particular polynomial
7
3
1
p( x )  x 4  3x 3  x 2  x 
4
4
2
(3) Now try your MATLAB program on the polynomial
12
p( x )   ( x  2  j )
j 0
Hint: Use Maple to rewrite this polynomial in standard form
n
p( x )   a j x j
j 0
M. Dumbser
17 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Efficient Evaluation of Polynomials
The standard format of defining a polynomial is computationally not the most efficient way to
evaluate a polynomial:
p( x )  a0  a1  x  a2  x 2  a3  x 3  a4  x 4 ...
p( x)  a0  a1  x  a2  x  x  a3  x  x  x  a4  x  x  x  x...
In this example, we have 4 additions and 10 multiplications
A computationally more efficient but mathematically fully equivalent way of writing the polynomial
p(x) is the so called Horner scheme, which uses the distributive law:
p( x)  a0  x  a1  x  a2  x  a3  x  a4 ...
Now, we still have 4 additions but only 4 multiplications!
M. Dumbser
18 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Globally Convergent Newton Methods
Exercise 4
Now try your original MATLAB script Newton.m to solve the nonlinear equation
f ( x )  arctan( x  1)  0
1)
Use x0 = 2 as initial guess.
2)
Use x0 = 2.5 as initial guess. What do you realize?
M. Dumbser
19 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Newton‘s Method for Systems of Equations
The extension of Newton‘s method to systems of nonlinear algebraic equations with several unknown
variables is straightforward.
The problem is written in standard form in matrix-vector notation as follows:
 
f ( x)  0
Linearization (multi-dimensional Taylor series up to first order) leads to:

 
f  
f ( xi )    x  xi   0
x
  
x : x  xi

f
f
D  Dij :   i
x x j
The matrix D is called the Jacobian matrix of the vector function f with respect to ist vector
argument x and is the multi-dimensional equivalent to the first derivative of a scalar function.
 

1
x   D f ( xi )



xi 1  xi  x
(1)
(2)
M. Dumbser
20 / 23
This corresponds to solving the linear system
 

Dx   f ( xi )
If ||f(xi)|| < tolerance then STOP (in Maple/MATLAB: break) or
If || x || < tolerance then STOP
 
f ( x) 
f12  f 22  ...  f n2
i.e. the vector norm of f is the so-called function residual.
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 6
Write a MATLAB script NewtonSys.m that solves the nonlinear system of equations
f1 ( x1, x2)  exp( x1)  sin( x2)  x12  x22  0
f 2 ( x1, x2)  cos( x1) sin( x2)  x1  x2  0
using the Newton method for systems. The user should have the possibility to choose
the initial guess for the vector of unknowns, the maximum number of iterations as
well as the tolerance. The vector function and its Jacobian matrix should be
implemented in two separate .m files, called FuncSys.m and DFuncSys.m.
M. Dumbser
21 / 23
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Globally Convergent Newton‘s Method for Systems
The extension of the globally convergent Newton‘s method to systems of nonlinear algebraic
equations with several unknown variables is done as in the scalar case:
DO i=0,1,…
 1
 f   

x     f ( xi )
 x 
0  1
WHILE
k 0
 
 

f ( xi   k xi )  f ( xi )
 k 1     k
k  k 1
END



xi 1  xi   k x
 
UNTIL f ( xi )  tol
M. Dumbser
22 / 23
0  1
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 7
We want to solve the following nonlinear system of equations:
f ( x, y )  arctan( x  1)  arctan( y )  0
g ( x, y )  exp(  x 2 ) 
1
y0
3
(1) Use the classical Newton method for systems starting from [0,0]
(2) Use the classical Newton method for systems starting from [1000,1000]
(3) Can you find a starting point from which the classical Newton method converges?
(4) Construct a globally convergent Newton method for systems and start again
from the points [0,0] and [1000,1000]. Does the global method even converge
from the starting point [1e6,1e6] ?
M. Dumbser
23 / 23