Transcript Day 4 Slides

Scientific Computing
Algorithm Convergence and
Root Finding Methods
Algorithm Convergence
• Definition: We say that a numerical algorithm to
solve some problem is convergent if the numerical
solution generated by the algorithm approaches the
actual solution as the number of steps in the
algorithm increases.
• Definition: Stability of an algorithm refers to the
ability of a numerical algorithm to reach the correct
solution, even if small errors are introduced during
the execution of the algorithm.
Algorithm Stability
• Stable algorithm: Small changes in the initial
data produce small changes in the final result
(Also, called Insensitive or Well-Conditioned)
• Unstable or conditionally stable algorithm:
Small changes in the initial data may produce
large changes in the final result
(Also called Sensitive or Ill-Conditioned)
Algorithm Stability
• Condition Number: The condition number of
a problem is the ratio of the relative change in
a solution to the relative change of the input
^
cond no 
y y
y
^
x x
x
y / y

x / x
• Want cond no to be small!
Algorithm Stability
• Example: Consider
p(x) = x5 – 10x4 + 40x3 – 80x2 + 80x – 32.
Note that p(x) = (x –2)5
Does Matlab roots command use a stable algorithm?
Algorithm Stability
• Example: Use the fzero Matlab function.
Does Matlab
fzero command
use a stable
algorithm?
Algorithm Error Analysis
• To calculate convergence, we usually try to
estimate the error in the numerical solution
compared to the actual solution. We try to
find out how fast this error is decreasing (or
increasing!)
Algorithm Error Disaster!!
• Explosion of the Ariane 5
• On June 4, 1996 an unmanned Ariane 5 rocket launched by
the European Space Agency exploded just forty seconds after
lift-off. The cause of the failure was a software error in the
inertial guidance system.
• A 64 bit floating point number relating to the horizontal
velocity of the rocket with respect to the platform was
converted to a 16 bit signed integer. The number was larger
than 32,768, the largest integer storeable in a 16 bit signed
integer, and thus the conversion failed.
(source http://ta.twi.tudelft.nl/users/vuik/wi211/disasters.html)
Algorithm Rate of Convergence
• To get an idea of how fast an algorithm
convergences, we measure the rate at which the
sequence of approximations generated by the
algorithm converge to some value.
Definition:
Algorithm Rate of Convergence
• Alternatively, we can consider how fast the error
goes to zero from one iteration of the algorithm to
the next.
• Definition: Let εn be the error in the nth iteration of
the algorithm. That is, εn = |αn - α |.
• Definition: We say that the algorithm has linear
convergence if |n+1| = C |n|
• Definition: We say that the algorithm has quadratic
convergence if |n+1| = C |n|2
• Which convergence type is faster?
Finding Roots (Zeroes) of Functions
• Given some function f(x) , find location x
where f(x)=0. This is called a root (or zero) of
f(x)
• For success we will need:
– Starting position x0, hopefully close to solution x
– Well-behaved function f(x)
Example
Problem: Determine the depth of an object in water without
submerging it.
x = depth of sphere
in water
1
1-x
1
1-x
x
r
r
Example
Problem: Determine the depth of an object in water without
submerging it.
Sample material: Cork –> ρ/ρw = 0.25 (Specific Gravity)
Solution: We found that the solution required finding the zero
of
3
2
x

3
x

1

0
Finding Roots (Zeroes) of Functions
• Notes:
– There may be multiple roots of f(x). That is why
we need to specify an initial guess of x0 .
– If f(x) is not continuous, we may miss a root no
matter what root-finding algorithm we try.
– The roots may be real numbers or complex
numbers. In this course we will consider only
functions with real roots.
Finding Roots (Zeroes) of Functions
• What else can go wrong?
Tangent point:
very difficult
to find
Singularity:
brackets don’t
surround root
Pathological case:
infinite number of
roots – e.g. sin(1/x)
Bisection method
• Simplest Root finding method. Based on
Intermediate Value Theorem: If f(x) is
continuous on [a, b] then for any y such that y
is between f(a) and f(b) there is a c in [a, b]
such that f(c) = y.
Bisection method
• Given points a and b that bracket a root, find
x = ½ (a+ b)
and evaluate f(x)
• If f(x) and f(b) have the same sign (+ or -) then
b  x else a  x
• Stop when a and b are “close enough”
• If function is continuous, this will succeed
in finding some root.
Bisection Matlab Function
function v = bisection(f, a, b, eps)
% Function to find root of (continuous) f(x) within [a,b]
% Assumes f(a) and f(b) bracket a root
k = 0;
while abs(b-a) > eps*abs(b)
x = (a + b)/2;
if sign(f(x)) == sign(f(b))
b = x;
else
a = x;
end
k = k + 1; root = x;
v = [root k];
end
Bisection Matlab Function
f = inline('x^3-3*x^2+1')
bisection(f, 1, 2, 0.01)
ans =
0.6523 7.0000
bisection(f, 1, 2, 0.005)
ans =
0.6543 8.0000
bisection(f, 1, 2, 0.0025)
ans =
0.6533 9.0000
bisection(f, 1, 2, 0.00125)
ans =
0.6528 10.0000
Is there a pattern?
Bisection Convergence
• Convergence rate:
– Error bounded by size of [a… b] interval
– Interval shrinks in half at each iteration
– Therefore, error cut in half at next iteration:
|n+1| = ½ |n|
– This is linear convergence
Root Finding Termination
• Let pn be the approximation generated by the
nth iteration of a root-finding algorithm. How
do we know when to stop?
• Bisection – uses a version of the second
condition.
Practice
• Class Exercise: Determine the number of iterations
needed to achieve an approximation to the solution
of x3 – 2x -1 = 0 lying within the interval [1.5, 2.0] and
accurate to 10-2 .
– Using the Bisection Method
• Estimate using theory, then try using Matlab