Introduction to Scientific
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter 8 Nonlinear Equations
and Optimization
? Finding Roots
? Minimizing a function of one variable
? Minimizing multivariate functions
? Solving systems of nonlinear equations
We consider several types of nonlinear problems,They differ in
whether or not the solution sought is a vector or a scalar and
whether or not the goal is to produce a root or a minimizer of
some given function,The presentation is organized around a
family of,orbit” problems,Suppose the vector-valued functions
??
?
??
??
??
?
??
??
)(
)()(
)(
)()(
2
2
2
1
1
1 ty
txtP
ty
txtP
Specify the location at time t of a pair of planets that go around
the Sun,Assume that the orbits are elliptical and that the sun is
situated at (0,0).
Question 1,At what times are the planets and the sun collinear?
If f(t) is the the sine of the angle between P1 and P2,then this
problem is equivalent to finding a zero of f(t),We focus on the
bisection and Newton methods,and the Matlab zero-finder,
fzero
Question 2,How close do the two planets get for a period of time?
If f(t) is distance from P1 to P2,then this is a single-variable
minimization problem,We develop the method of golden section
search and discuss the Matlab minmizer fmin.
Question 3,How close do the two orbits get? The method of
steepest descent and Matlab multivariable minimizer fmins are
designed to solve problems of this variety.
Question 4,Where (if at all) do the two orbits intersect? This is an
example of a multivariable root-finding problem:
The Newton framework for systems of nonlinear equations is
discussed,Related topics include the use of finite differences to
approximate the Jacobian and the Gauss-Newton method for the
nonlinear least squares problem.
0)()(),( 112221 ??? tPtPttF
22211
m i n
,)()(21 tPtPtt ?
Finding Roots
06210)( 53 ???? ?? tt eetf
quintic polynomial 072 245 ????? xxxx
Algorithms in this area are iterative and proceed by producing a
sequence of numbers that converge to a root of interest.
Where do we start the iteration?
Does the iteration converge and how fast?
How do we know when to terminate the iteration?
The square root problem
Suppose A is a positive real number and that we want to
compute its square root,Geometrically,this task is equivalent
to constructing a square whose area is A.
??
?
?
??
?
?
???
c
c x
A
xx
2
1
xc 6000/xc
------------------------------------------- ------- ------
60.00000000000000 100.00000000000000
80.00000000000000 75.00000000000000
77.50000000000000 77.41935483870968
77.45967741935485 77.45965642894325
77.45966692414905 77.45966692414763
77.45966692414834 77.45966692414834
77.45966692414834 77.45966692414834
77.45966692414834
To improve an approximate square root
141,4* ??? mmA e
]1,25[.,2* ?? mmA e
a good initial guess for the reduced range problem can be obtained
by linear interpolation with
3/)21()( mmL ??
05.0|)(| ?? mmL
0 0.5 1 1.5
0
0.2
0.4
0.6
0.8
1
1.2
m
sqrt(m)
2
2
1
2
1
?
?
?
?
?
?
?
? ?
????
?
?
??
?
?
????
c
c
c
c
x
mx
m
x
m
xmx
mxemLx kk ??? ),(0 21 2 1 k
k
k exe ??
It can be shown that the iterates xk are always in the interval
[.5,1],and so
21 kk ee ??
161608142234 )05(.????? eeeee
function x = MySqrt(A)
if A==0,x = 0; else
TwoPower = 1; m = A;
while m>=1,m = m/4; TwoPower = 2*TwoPower;
end
while m <,25
m = m*4; TwoPower = TwoPower/2; end
x = (1+2*m)/3;
for k=1:4,x = (x + (m/x))/2; end
x = x*TwoPower;
end
0 2 4 6 8 10 12 14 160
0.5
1
1.5
2
2.5
3 x 10
-16 Relative Error in the Function MySqrt(A)
A
Bisection
If a continuous function f(x) changes sign on an interval,then
it must have at least one root in the interval.
while abs(a-b) > delta
if f(a)*f((a+b)/2)<=0
b = (a+b)/2
else
a = (a+b)/2
end
end
root = (a+b)/2;
The while-loop may never
terminate if delta is smaller
than the spacing of the
floating point numbers
between a and b.
The second flaw is that it
requires two function
evaluations per iteration.
function noend(a,b,delta)
i=1;
while abs(a-b)>delta
if a<=b,a=(a+b)/2; else,b=(a+b)/2;
end
i=i+1
if abs(a-b)<eps*max(abs(a),abs(b))
a
pause
end
end
noend(1000.0001,1001.0001,eps)
noend(10000000.0001,10000001.0001,eps)
function root = Bisection(fname,a,b,delta)
fa = feval(fname,a); fb = feval(fname,b);
if fa*fb > 0,disp('Initial interval is not bracketing.')
return;end
if nargin==3,delta = 0;end
while abs(a-b) > delta+eps*max(abs(a),abs(b))
mid = (a+b)/2; fmid = feval(fname,mid);
if fa*fmid<=0,b = mid; fb = fmid;
else,a = mid; fa = fmid; end
end
root = (a+b)/2;
kkk
baba
2
|||| 00 ???
1
00
* 2
||||,
2
)(
?
?????
kk
kk
k
baxxbax
Linear convergence
)1,0[|,||| **1 ????? cxxcxx kk
myshowBisection.m
Showbisect.m
Notice that a new digit of pi is acquired
every three or so iterations,This kind of
uniform acquisition of significant digits is
the hallmark of methods that converge
linearly
1 1, 5 2 2, 5 3 3, 5 4 4, 5 5 5, 5 6
-2
0
2
4
6
8
10
12
14
x c = 5, 2 6 8 4 3 3 f c = 2, 8 5 7 e + 0 0 0
The Newton Method Idea
)(')()()( cccc xfxxxfxL ???
)('
)(
c
c
c xf
xfxx ??
?
k xk xk-pi
-------------------------------------------------------------------
0 3.1415926535897931 0.0000000000000000
1 3.7963140465723391 0.6547213929825460
2 3.2594354361754703 0.1178427825856772
3 3.1451315542075164 0.0035389006177233
4 3.1415957863900608 0.0000031328002676
5 3.1415926535922467 0.0000000000024536
6 3.1415926535897931 0.0000000000000000
1)4/t a n ()( ?? xxf
Mytestnewton.m
quadratical convergence
,|||| 2**1 xxcxx kk ????
No guarantee that x+ is closer to a root than xc
A small value for f’(x *) works against rapid convergence
0,10,10,10,10,2,)( 1284022 ??????? ataxxf
||||,3 9 1 7.1||
,
1
1
)('),a r c t a n ()(
2
cc xxxif
x
xfxxf
??
?
??
?
divergence
q u a d r a t ic is ec o n v e r g e n c r a t e T h e,i n v a l u es t a r t i n g
a n yf o r ec o n v e r g e n c ngg u a r a n t e e it h e r e b y,t od i s t a n c e
t h eh a l v e s t h a n m o r e s t e pN e w t o n a,a n y f o r T h u s
||
2
1
||
2
||
a n d,)(')/( t h e n if
/
.,a l lf o r |y-x||)(')('|
,a l lf o r |)('|
a t hs u c h e x i s t a n d
c o n s t a n t s p o s i t i v e a n d 0,0)( w h e r e,],[
i n t e r v a la n o n d e f i n e d a r e )(' a n d )( s u p p o s e 8 T h e o r e m
*
*
2
**
***
I
x
Ix
xxxxxx
IxfxfxxIx
Iyxyfxf
Ixxf
xfxxI
xfxf
c
cc
cccc
?
?????
????
?
???
??
?????
?
?
?
?
???
?
?
??
???
Globalizing the Newton Iteration
By combining the Newton and bisection ideas in a way that
captures the best features of each framework.
Holt condition
The length of the current bracketing interval is less that tolx,This
ensures that the computed root is within tolx of a true root,It does
not guarantee that f is small.
The absolue value of f(xc) is less than tolf,This does not mean that
xc is close to a true root.
The number of f-evaluations exceeds a specified positive integer
nEvalsMax,This means that neither of the preceding two conditions
is satisfied.
)(,0 0 0 0 0 0 0 0 0 1.1,)1()( 1,???? cc xfxxxf
)(,1.1,)1()( 10 ???? cc xfxxxf
step a x b
------------------------------------------------------------------------------------------------
<start> -10.995574287564276 -10.995574287564276 47.223889803846895
Bisection -10.995574287564276 18.114157758141310 18.114157758141310
Bisection -10.995574287564276 3.559291735288517 3.559291735288517
Newton 3.115476144648328 3.115476144648328 3.559291735288517
Newton 3.115476144648328 3.141598592990409 3.141598592990409
Newton 3.141592653589793 3.141592653589793 3.141598592990409
Avoiding Derivatives—Secant method
,)()()('
?
?
?
??
xx
xfxfxf
c
c
c
6.12 51,|||| **1 ??????? rxxcxx rkk
0 1.0000000000000000 -2.1415926535897931
1 3.5593092641513646 0.4177166105615715
2 2.9470453609190379 -0.1945472926707552
3 3.1208629615047401 -0.0207296920850530
4 3.1425827065470964 0.0009900529573033
5 3.1415875311545123 -0.0000051224352808
6 3.1415926523218185 -0.0000000012679746
7 3.1415926535897949 0.0000000000000018
8 3.1415926535897931 0.0000000000000000
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8-6
-5
-4
-3
-2
-1
0
1
2
3
4
The Matlab fzero function
Function y=decay(t)
y=10*exp(-3*t)+2*exp(-2*t)-6;
Root=fzero(?decay?,1);
close all
global fig;global mov;
global k;k=0;
fig=figure(1);
axis([0,2,-7,4]);
hold on
set(gca,'xlim',[0 2],'ylim',[-7 4],'NextPlot','add','Visible','on')
mov = avifile('findzero1.avi','fps',5)
root=fzero('decay',1);
mov=close(mov);
function y=decay(t)
global k; global fig;global mov
y=10*exp(-3*t)+2*exp(-2*t)-6;
figure(fig);axis([0,2,-7,4]);
plot(t,y,'*');k=k+1;
text(t,y,sprintf('%3d',k));
F = getframe(fig);
for t=1:5
mov = addframe(mov,F);
end
An application
)97.87/2s i n (6 7 4 1.56)(
)97.87/2c o s (9 1 1 7.579 0 8 4.11)(
tty
ttx
M
M
?
?
?
???
)25.3 6 5/2s i n (5 8 3 2.1 4 9)(
)25.3 6 5/2c o s (6 0 4 1.1 4 94 9 8 7.2)(
tty
ttx
E
E
?
?
?
???
Mercury
Earth
To an observer on E,M is in conjunction if it is located
on the Sun-to-Earth line segment.
0 50 100 150 200 250 300 350 400-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Sine of the Mercury-Sun-Earth Angle
Next ten conjunctions,Spacing Estimate = 115.000
Conjunction Time Spacing
1 112.476 112.476
2 234.682 122.206
3 348.554 113.872
4 459.986 111.433
5 581.491 121.505
6 697.052 115.561
7 807.815 110.762
8 928.020 120.206
9 1045.440 117.420
10 1155.908 110.467
Minimizing a Function of a single Variable
}1 0 0 00:{),(m i n ???? ttStfSt
22 ))()(())()(()( tytytxtxtf MEME ????
Optimization problem
Objective function
0 100 200 300 400 500 600 700 800 900 1000
80
100
120
140
160
180
200
220 The Function distmercearth
Enter New Left Endpoint,(Click off x-range to terminate.)
83.3061
Graphical Search
plot(x,y,[L R],[ymin ymin])
[x1,y1] = ginput(1);
function [L,R] = GraphSearch(fname,a,b,Save,nfevals)
870 880 890 900 910 920 930 940 950 960 970
80
100
120
140
160
180
200
220 The Function distmercearth
Enter New Left Endpoint,(Click off x-range to terminate.)
82.6812
Strictly
monotone
decreasin
g
Strictly
monotone
increasing
unimodal
c da b
Golden Section Section
c da b
))(1(
)(
abrad
abrac
????
???
Golden
section search
0<r<1/2,If f(c) ≤f(d),then from unimodality it follows that f
is increasing in [d,b] and so t* in [a,d],otherwise t* is in [c,d]
After this step the length of the search interval is reduced by a
factor of (1-r).
))(1()(
)())(1(
abraabrc
abraadra
??????
??????
2
53 ??r
6 1 8.11
*
??? rrr
* is golden ration
Golden Section Search,(a c d b) Trace
0
Step
1
2
3
4
5
golden('distmercearth',900,950)
927.1243
Termination criteria in Golden
while (d-c) > sqrt(eps)*max(abs(c),abs(d))
)(''2)()( *
2
** xfxfxf
?? ???
)(''2)( *
2
* xfxf
?? ??
Unimodality is to Golden as the requirement f(a)f(b) ≤0
The Matlab fmin function
The Matlab minimizer fmin uses a combination of golden
section search and a parabolic fitting method that relies only
upon f evalution.
tmin = fmin(?DistMercEarth?,900,950)
Warning,FMIN is obsolete and has been replaced by FMINBND.
FMIN now calls FMINBND which uses the following syntax:
[X,FVAL,EXITFLAG,OUTPUT] =
FMINBND(FUN,x1,x2,OPTIONS,P1,P2,...)
Use OPTIMSET to define optimization options,or type
'edit fmin' to view the code used here,FMIN will be
removed in the future; please use FMINBND with the new syntax.
Minimizing Multivariate
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
)s in (
)c o s (
22
)c o s ()s in (
)s in ()c o s (
)(
)(
11
1111
11
11
1
1
tAP
t
APAP
ty
tx
??
??
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
)s in (
)c o s (
22
)c o s ()s in (
)s in ()c o s (
)(
)(
22
2222
221
22
2
2
tAP
t
APAP
ty
tx
??
??
? ?222112221121 ))()(())()((21),( tytytxtxtts e p ????
-10 -8 -6 -4 -2 0 2 4-4
-2
0
2
4
6
8
min Sep = 0.0785
function pLoc = Orbit(t,planet,lineStyle)
c = cos(t); s = sin(t);
x0 = ((planet.P-planet.A)/2) + ((planet.P+planet.A)/2)*c;
y0 = sqrt(planet.A*planet.P)*s;
cphi = cos(planet.phi); sphi = sin(planet.phi);
pLoc = struct('x',cphi*x0 + sphi*y0,'y',-sphi*x0 + cphi*y0);
if nargin==3
plot(pLoc.x,pLoc.y,lineStyle)
end
The Gradient
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
n
c
c
c
t
tf
t
tf
tf
)(
)(
)(
1
?
?
?
?
?
?
?
????
???
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
)()]()([)()]()([
)()]()([)()]()([
),(
),(
)(
222211222211
112211112211
2
21
1
21
tytytytxtxtx
tytytytxtxtx
t
tts e p
t
tts e p
ts e p
??
??
0 1 2 3 4 5 60
1
2
3
4
5
6
Enter Contour Labels (Optional),To quit,strike,
t1
t2
57.5
23.1
5.83
28.8
5.83
46
The Method of Steepest Descent
)()()()( ccTcccccc Ostftfstf ??? ?????
)( cc tfs ???
ccc
c
cccc
cccc
stt
tf
stf
gstfg
?
??
???
??
?????
?
S e t
)( t h a n le s sly s u f f ic ie n t is
)( so c h o o s ee ly I n e x p e n s iv
s e t a n d )(g r a d ie n t t h ec o m p u t e
0 1 2 3 4 5 6
0
1
2
3
4
5
6
tmin = ( 5.241,3.926) norm(gmin)= 3.4818e-001
t1
t2
57.5
23.1
5.83
28.8
5.83
46
0 5 10 15 20 25 30
0
5
10
15
20
Value of Sep
Iteration
0 5 10 15 20 25 30
10-1
100
101
102
Value of norm(gSep).
Iteration
-10 -8 -6 -4 -2 0 2 4-4
-2
0
2
4
6
8
min Sep = 0.1277
0 1 2 3 4 5 6
0
1
2
3
4
5
6
t m i n = ( - 0, 9 7 6,4, 0 3 8 ) n o r m ( g m i n ) = 3, 5 7 7 0 e - 0 0 1
t1
t2
5 7, 5
1 7, 3
5, 8 3
5 1, 8
5, 8 3
2 3, 1
1 1, 6
0 5 10 15 20 25 300
10
20
30
40
50
60
Value of Sep
Iteration
0 5 10 15 20 25 3010
-1
100
101
102
Value of norm(gSep).
Iteration
-10 -8 -6 -4 -2 0 2 4-4
-2
0
2
4
6
8
min Sep = 0.1005
The Matlab Fmins Function
The fmins function can be used to minimize a multivariate
function,It uses a simplex search method and does not require
gradient function implementation.
Warning,FMINS is obsolete and has been replaced by
FMINSEARCH.
FMINS now calls FMINSEARCH which uses the following
syntax:
[X,FVAL,EXITFLAG,OUTPUT] =
FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...)
Use OPTIMSET to define optimization options,or type
'edit fmins' to view the code used here,FMINS will be
removed in the future; please use FMINSEARCH with the new
syntax.
t* = ( 5.501623,4.467533 ),Steps = 103,norm(gradient) = 1.915e-006
Systems of Nonlinear equations
-20 -15 -10 -5 0 5-12
-10
-8
-6
-4
-2
0
2
4
6
8
Initial Guess
0
)()(
)()(
)(
1122
112 ?
??
?
??
?
?
?
?
tyty
txtx
ts e p V
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
)s in (
)c o s (
22
)c o s ()s in (
)s in ()c o s (
)(
)(
tAP
t
APAP
ty
tx
ii
iiii
ii
ii
i
i
??
??
ccccc stJtFstF )()()( ???
??
?
??
?
?
?
???
)()(
)()(
)()(
2211
2211
tyty
txtx
tFtJ cc ??
??
The Newton Idea
0)()()( ???? ccccc stJtFstF
)()( ccc tFstJ ??
Nonlinear Least Squares
?
?
?
?
?
?
?
? ???
??
?
?
?
?
?
)s in (
)c o s (
22
)(
)(
tAP
t
PAAP
ty
tx
))c o s (1())c o s (1(
2)(
??? ???? AP
APr
))c o s (1())c o s (1(
2),(
ii
ii AP
APrPAF
?? ?????
?
?
?
m
i
i PAFPAO r b
1
2),(
2
1),(minimize
?
?
?
?
?
?
?
?
?
?
?
?
?
)(
)(
)(
)(
2
1
pF
pF
pF
pF
m
?
?
?
?
m
i
i pFp
1
2)(
2
1)(?
To seek a vector of unknown parameters p so that the sum
of squares is as small as possible
)()(2 ccc psp ?? ????
A full Newton method involving the step
entails derivative evaluation since the gradient is given by
?
?
????
??
1
22 )()()()()(
)()()(
i
cicic
T
cc
c
T
cc
pFpFpJpJp
pFpJp
?
?
Gauss-Newton framework
cc
c
T
ccc
T
c
cccc
spp
FJsJJ
pJJpFF
???
???
???
? U p d a t e
)( S o l v e
)(J a c o b i a n t h ea n d )( E v a l u a t e
2)()(m i n ccRs pFspJn ??
Note that sc is the solution to the normal equations for the least
squares problem
-5 -4 -3 -2 -1 0 1
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter 8 Nonlinear Equations
and Optimization
? Finding Roots
? Minimizing a function of one variable
? Minimizing multivariate functions
? Solving systems of nonlinear equations
We consider several types of nonlinear problems,They differ in
whether or not the solution sought is a vector or a scalar and
whether or not the goal is to produce a root or a minimizer of
some given function,The presentation is organized around a
family of,orbit” problems,Suppose the vector-valued functions
??
?
??
??
??
?
??
??
)(
)()(
)(
)()(
2
2
2
1
1
1 ty
txtP
ty
txtP
Specify the location at time t of a pair of planets that go around
the Sun,Assume that the orbits are elliptical and that the sun is
situated at (0,0).
Question 1,At what times are the planets and the sun collinear?
If f(t) is the the sine of the angle between P1 and P2,then this
problem is equivalent to finding a zero of f(t),We focus on the
bisection and Newton methods,and the Matlab zero-finder,
fzero
Question 2,How close do the two planets get for a period of time?
If f(t) is distance from P1 to P2,then this is a single-variable
minimization problem,We develop the method of golden section
search and discuss the Matlab minmizer fmin.
Question 3,How close do the two orbits get? The method of
steepest descent and Matlab multivariable minimizer fmins are
designed to solve problems of this variety.
Question 4,Where (if at all) do the two orbits intersect? This is an
example of a multivariable root-finding problem:
The Newton framework for systems of nonlinear equations is
discussed,Related topics include the use of finite differences to
approximate the Jacobian and the Gauss-Newton method for the
nonlinear least squares problem.
0)()(),( 112221 ??? tPtPttF
22211
m i n
,)()(21 tPtPtt ?
Finding Roots
06210)( 53 ???? ?? tt eetf
quintic polynomial 072 245 ????? xxxx
Algorithms in this area are iterative and proceed by producing a
sequence of numbers that converge to a root of interest.
Where do we start the iteration?
Does the iteration converge and how fast?
How do we know when to terminate the iteration?
The square root problem
Suppose A is a positive real number and that we want to
compute its square root,Geometrically,this task is equivalent
to constructing a square whose area is A.
??
?
?
??
?
?
???
c
c x
A
xx
2
1
xc 6000/xc
------------------------------------------- ------- ------
60.00000000000000 100.00000000000000
80.00000000000000 75.00000000000000
77.50000000000000 77.41935483870968
77.45967741935485 77.45965642894325
77.45966692414905 77.45966692414763
77.45966692414834 77.45966692414834
77.45966692414834 77.45966692414834
77.45966692414834
To improve an approximate square root
141,4* ??? mmA e
]1,25[.,2* ?? mmA e
a good initial guess for the reduced range problem can be obtained
by linear interpolation with
3/)21()( mmL ??
05.0|)(| ?? mmL
0 0.5 1 1.5
0
0.2
0.4
0.6
0.8
1
1.2
m
sqrt(m)
2
2
1
2
1
?
?
?
?
?
?
?
? ?
????
?
?
??
?
?
????
c
c
c
c
x
mx
m
x
m
xmx
mxemLx kk ??? ),(0 21 2 1 k
k
k exe ??
It can be shown that the iterates xk are always in the interval
[.5,1],and so
21 kk ee ??
161608142234 )05(.????? eeeee
function x = MySqrt(A)
if A==0,x = 0; else
TwoPower = 1; m = A;
while m>=1,m = m/4; TwoPower = 2*TwoPower;
end
while m <,25
m = m*4; TwoPower = TwoPower/2; end
x = (1+2*m)/3;
for k=1:4,x = (x + (m/x))/2; end
x = x*TwoPower;
end
0 2 4 6 8 10 12 14 160
0.5
1
1.5
2
2.5
3 x 10
-16 Relative Error in the Function MySqrt(A)
A
Bisection
If a continuous function f(x) changes sign on an interval,then
it must have at least one root in the interval.
while abs(a-b) > delta
if f(a)*f((a+b)/2)<=0
b = (a+b)/2
else
a = (a+b)/2
end
end
root = (a+b)/2;
The while-loop may never
terminate if delta is smaller
than the spacing of the
floating point numbers
between a and b.
The second flaw is that it
requires two function
evaluations per iteration.
function noend(a,b,delta)
i=1;
while abs(a-b)>delta
if a<=b,a=(a+b)/2; else,b=(a+b)/2;
end
i=i+1
if abs(a-b)<eps*max(abs(a),abs(b))
a
pause
end
end
noend(1000.0001,1001.0001,eps)
noend(10000000.0001,10000001.0001,eps)
function root = Bisection(fname,a,b,delta)
fa = feval(fname,a); fb = feval(fname,b);
if fa*fb > 0,disp('Initial interval is not bracketing.')
return;end
if nargin==3,delta = 0;end
while abs(a-b) > delta+eps*max(abs(a),abs(b))
mid = (a+b)/2; fmid = feval(fname,mid);
if fa*fmid<=0,b = mid; fb = fmid;
else,a = mid; fa = fmid; end
end
root = (a+b)/2;
kkk
baba
2
|||| 00 ???
1
00
* 2
||||,
2
)(
?
?????
kk
kk
k
baxxbax
Linear convergence
)1,0[|,||| **1 ????? cxxcxx kk
myshowBisection.m
Showbisect.m
Notice that a new digit of pi is acquired
every three or so iterations,This kind of
uniform acquisition of significant digits is
the hallmark of methods that converge
linearly
1 1, 5 2 2, 5 3 3, 5 4 4, 5 5 5, 5 6
-2
0
2
4
6
8
10
12
14
x c = 5, 2 6 8 4 3 3 f c = 2, 8 5 7 e + 0 0 0
The Newton Method Idea
)(')()()( cccc xfxxxfxL ???
)('
)(
c
c
c xf
xfxx ??
?
k xk xk-pi
-------------------------------------------------------------------
0 3.1415926535897931 0.0000000000000000
1 3.7963140465723391 0.6547213929825460
2 3.2594354361754703 0.1178427825856772
3 3.1451315542075164 0.0035389006177233
4 3.1415957863900608 0.0000031328002676
5 3.1415926535922467 0.0000000000024536
6 3.1415926535897931 0.0000000000000000
1)4/t a n ()( ?? xxf
Mytestnewton.m
quadratical convergence
,|||| 2**1 xxcxx kk ????
No guarantee that x+ is closer to a root than xc
A small value for f’(x *) works against rapid convergence
0,10,10,10,10,2,)( 1284022 ??????? ataxxf
||||,3 9 1 7.1||
,
1
1
)('),a r c t a n ()(
2
cc xxxif
x
xfxxf
??
?
??
?
divergence
q u a d r a t ic is ec o n v e r g e n c r a t e T h e,i n v a l u es t a r t i n g
a n yf o r ec o n v e r g e n c ngg u a r a n t e e it h e r e b y,t od i s t a n c e
t h eh a l v e s t h a n m o r e s t e pN e w t o n a,a n y f o r T h u s
||
2
1
||
2
||
a n d,)(')/( t h e n if
/
.,a l lf o r |y-x||)(')('|
,a l lf o r |)('|
a t hs u c h e x i s t a n d
c o n s t a n t s p o s i t i v e a n d 0,0)( w h e r e,],[
i n t e r v a la n o n d e f i n e d a r e )(' a n d )( s u p p o s e 8 T h e o r e m
*
*
2
**
***
I
x
Ix
xxxxxx
IxfxfxxIx
Iyxyfxf
Ixxf
xfxxI
xfxf
c
cc
cccc
?
?????
????
?
???
??
?????
?
?
?
?
???
?
?
??
???
Globalizing the Newton Iteration
By combining the Newton and bisection ideas in a way that
captures the best features of each framework.
Holt condition
The length of the current bracketing interval is less that tolx,This
ensures that the computed root is within tolx of a true root,It does
not guarantee that f is small.
The absolue value of f(xc) is less than tolf,This does not mean that
xc is close to a true root.
The number of f-evaluations exceeds a specified positive integer
nEvalsMax,This means that neither of the preceding two conditions
is satisfied.
)(,0 0 0 0 0 0 0 0 0 1.1,)1()( 1,???? cc xfxxxf
)(,1.1,)1()( 10 ???? cc xfxxxf
step a x b
------------------------------------------------------------------------------------------------
<start> -10.995574287564276 -10.995574287564276 47.223889803846895
Bisection -10.995574287564276 18.114157758141310 18.114157758141310
Bisection -10.995574287564276 3.559291735288517 3.559291735288517
Newton 3.115476144648328 3.115476144648328 3.559291735288517
Newton 3.115476144648328 3.141598592990409 3.141598592990409
Newton 3.141592653589793 3.141592653589793 3.141598592990409
Avoiding Derivatives—Secant method
,)()()('
?
?
?
??
xx
xfxfxf
c
c
c
6.12 51,|||| **1 ??????? rxxcxx rkk
0 1.0000000000000000 -2.1415926535897931
1 3.5593092641513646 0.4177166105615715
2 2.9470453609190379 -0.1945472926707552
3 3.1208629615047401 -0.0207296920850530
4 3.1425827065470964 0.0009900529573033
5 3.1415875311545123 -0.0000051224352808
6 3.1415926523218185 -0.0000000012679746
7 3.1415926535897949 0.0000000000000018
8 3.1415926535897931 0.0000000000000000
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8-6
-5
-4
-3
-2
-1
0
1
2
3
4
The Matlab fzero function
Function y=decay(t)
y=10*exp(-3*t)+2*exp(-2*t)-6;
Root=fzero(?decay?,1);
close all
global fig;global mov;
global k;k=0;
fig=figure(1);
axis([0,2,-7,4]);
hold on
set(gca,'xlim',[0 2],'ylim',[-7 4],'NextPlot','add','Visible','on')
mov = avifile('findzero1.avi','fps',5)
root=fzero('decay',1);
mov=close(mov);
function y=decay(t)
global k; global fig;global mov
y=10*exp(-3*t)+2*exp(-2*t)-6;
figure(fig);axis([0,2,-7,4]);
plot(t,y,'*');k=k+1;
text(t,y,sprintf('%3d',k));
F = getframe(fig);
for t=1:5
mov = addframe(mov,F);
end
An application
)97.87/2s i n (6 7 4 1.56)(
)97.87/2c o s (9 1 1 7.579 0 8 4.11)(
tty
ttx
M
M
?
?
?
???
)25.3 6 5/2s i n (5 8 3 2.1 4 9)(
)25.3 6 5/2c o s (6 0 4 1.1 4 94 9 8 7.2)(
tty
ttx
E
E
?
?
?
???
Mercury
Earth
To an observer on E,M is in conjunction if it is located
on the Sun-to-Earth line segment.
0 50 100 150 200 250 300 350 400-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Sine of the Mercury-Sun-Earth Angle
Next ten conjunctions,Spacing Estimate = 115.000
Conjunction Time Spacing
1 112.476 112.476
2 234.682 122.206
3 348.554 113.872
4 459.986 111.433
5 581.491 121.505
6 697.052 115.561
7 807.815 110.762
8 928.020 120.206
9 1045.440 117.420
10 1155.908 110.467
Minimizing a Function of a single Variable
}1 0 0 00:{),(m i n ???? ttStfSt
22 ))()(())()(()( tytytxtxtf MEME ????
Optimization problem
Objective function
0 100 200 300 400 500 600 700 800 900 1000
80
100
120
140
160
180
200
220 The Function distmercearth
Enter New Left Endpoint,(Click off x-range to terminate.)
83.3061
Graphical Search
plot(x,y,[L R],[ymin ymin])
[x1,y1] = ginput(1);
function [L,R] = GraphSearch(fname,a,b,Save,nfevals)
870 880 890 900 910 920 930 940 950 960 970
80
100
120
140
160
180
200
220 The Function distmercearth
Enter New Left Endpoint,(Click off x-range to terminate.)
82.6812
Strictly
monotone
decreasin
g
Strictly
monotone
increasing
unimodal
c da b
Golden Section Section
c da b
))(1(
)(
abrad
abrac
????
???
Golden
section search
0<r<1/2,If f(c) ≤f(d),then from unimodality it follows that f
is increasing in [d,b] and so t* in [a,d],otherwise t* is in [c,d]
After this step the length of the search interval is reduced by a
factor of (1-r).
))(1()(
)())(1(
abraabrc
abraadra
??????
??????
2
53 ??r
6 1 8.11
*
??? rrr
* is golden ration
Golden Section Search,(a c d b) Trace
0
Step
1
2
3
4
5
golden('distmercearth',900,950)
927.1243
Termination criteria in Golden
while (d-c) > sqrt(eps)*max(abs(c),abs(d))
)(''2)()( *
2
** xfxfxf
?? ???
)(''2)( *
2
* xfxf
?? ??
Unimodality is to Golden as the requirement f(a)f(b) ≤0
The Matlab fmin function
The Matlab minimizer fmin uses a combination of golden
section search and a parabolic fitting method that relies only
upon f evalution.
tmin = fmin(?DistMercEarth?,900,950)
Warning,FMIN is obsolete and has been replaced by FMINBND.
FMIN now calls FMINBND which uses the following syntax:
[X,FVAL,EXITFLAG,OUTPUT] =
FMINBND(FUN,x1,x2,OPTIONS,P1,P2,...)
Use OPTIMSET to define optimization options,or type
'edit fmin' to view the code used here,FMIN will be
removed in the future; please use FMINBND with the new syntax.
Minimizing Multivariate
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
)s in (
)c o s (
22
)c o s ()s in (
)s in ()c o s (
)(
)(
11
1111
11
11
1
1
tAP
t
APAP
ty
tx
??
??
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
)s in (
)c o s (
22
)c o s ()s in (
)s in ()c o s (
)(
)(
22
2222
221
22
2
2
tAP
t
APAP
ty
tx
??
??
? ?222112221121 ))()(())()((21),( tytytxtxtts e p ????
-10 -8 -6 -4 -2 0 2 4-4
-2
0
2
4
6
8
min Sep = 0.0785
function pLoc = Orbit(t,planet,lineStyle)
c = cos(t); s = sin(t);
x0 = ((planet.P-planet.A)/2) + ((planet.P+planet.A)/2)*c;
y0 = sqrt(planet.A*planet.P)*s;
cphi = cos(planet.phi); sphi = sin(planet.phi);
pLoc = struct('x',cphi*x0 + sphi*y0,'y',-sphi*x0 + cphi*y0);
if nargin==3
plot(pLoc.x,pLoc.y,lineStyle)
end
The Gradient
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
n
c
c
c
t
tf
t
tf
tf
)(
)(
)(
1
?
?
?
?
?
?
?
????
???
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
)()]()([)()]()([
)()]()([)()]()([
),(
),(
)(
222211222211
112211112211
2
21
1
21
tytytytxtxtx
tytytytxtxtx
t
tts e p
t
tts e p
ts e p
??
??
0 1 2 3 4 5 60
1
2
3
4
5
6
Enter Contour Labels (Optional),To quit,strike,
t1
t2
57.5
23.1
5.83
28.8
5.83
46
The Method of Steepest Descent
)()()()( ccTcccccc Ostftfstf ??? ?????
)( cc tfs ???
ccc
c
cccc
cccc
stt
tf
stf
gstfg
?
??
???
??
?????
?
S e t
)( t h a n le s sly s u f f ic ie n t is
)( so c h o o s ee ly I n e x p e n s iv
s e t a n d )(g r a d ie n t t h ec o m p u t e
0 1 2 3 4 5 6
0
1
2
3
4
5
6
tmin = ( 5.241,3.926) norm(gmin)= 3.4818e-001
t1
t2
57.5
23.1
5.83
28.8
5.83
46
0 5 10 15 20 25 30
0
5
10
15
20
Value of Sep
Iteration
0 5 10 15 20 25 30
10-1
100
101
102
Value of norm(gSep).
Iteration
-10 -8 -6 -4 -2 0 2 4-4
-2
0
2
4
6
8
min Sep = 0.1277
0 1 2 3 4 5 6
0
1
2
3
4
5
6
t m i n = ( - 0, 9 7 6,4, 0 3 8 ) n o r m ( g m i n ) = 3, 5 7 7 0 e - 0 0 1
t1
t2
5 7, 5
1 7, 3
5, 8 3
5 1, 8
5, 8 3
2 3, 1
1 1, 6
0 5 10 15 20 25 300
10
20
30
40
50
60
Value of Sep
Iteration
0 5 10 15 20 25 3010
-1
100
101
102
Value of norm(gSep).
Iteration
-10 -8 -6 -4 -2 0 2 4-4
-2
0
2
4
6
8
min Sep = 0.1005
The Matlab Fmins Function
The fmins function can be used to minimize a multivariate
function,It uses a simplex search method and does not require
gradient function implementation.
Warning,FMINS is obsolete and has been replaced by
FMINSEARCH.
FMINS now calls FMINSEARCH which uses the following
syntax:
[X,FVAL,EXITFLAG,OUTPUT] =
FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...)
Use OPTIMSET to define optimization options,or type
'edit fmins' to view the code used here,FMINS will be
removed in the future; please use FMINSEARCH with the new
syntax.
t* = ( 5.501623,4.467533 ),Steps = 103,norm(gradient) = 1.915e-006
Systems of Nonlinear equations
-20 -15 -10 -5 0 5-12
-10
-8
-6
-4
-2
0
2
4
6
8
Initial Guess
0
)()(
)()(
)(
1122
112 ?
??
?
??
?
?
?
?
tyty
txtx
ts e p V
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
)s in (
)c o s (
22
)c o s ()s in (
)s in ()c o s (
)(
)(
tAP
t
APAP
ty
tx
ii
iiii
ii
ii
i
i
??
??
ccccc stJtFstF )()()( ???
??
?
??
?
?
?
???
)()(
)()(
)()(
2211
2211
tyty
txtx
tFtJ cc ??
??
The Newton Idea
0)()()( ???? ccccc stJtFstF
)()( ccc tFstJ ??
Nonlinear Least Squares
?
?
?
?
?
?
?
? ???
??
?
?
?
?
?
)s in (
)c o s (
22
)(
)(
tAP
t
PAAP
ty
tx
))c o s (1())c o s (1(
2)(
??? ???? AP
APr
))c o s (1())c o s (1(
2),(
ii
ii AP
APrPAF
?? ?????
?
?
?
m
i
i PAFPAO r b
1
2),(
2
1),(minimize
?
?
?
?
?
?
?
?
?
?
?
?
?
)(
)(
)(
)(
2
1
pF
pF
pF
pF
m
?
?
?
?
m
i
i pFp
1
2)(
2
1)(?
To seek a vector of unknown parameters p so that the sum
of squares is as small as possible
)()(2 ccc psp ?? ????
A full Newton method involving the step
entails derivative evaluation since the gradient is given by
?
?
????
??
1
22 )()()()()(
)()()(
i
cicic
T
cc
c
T
cc
pFpFpJpJp
pFpJp
?
?
Gauss-Newton framework
cc
c
T
ccc
T
c
cccc
spp
FJsJJ
pJJpFF
???
???
???
? U p d a t e
)( S o l v e
)(J a c o b i a n t h ea n d )( E v a l u a t e
2)()(m i n ccRs pFspJn ??
Note that sc is the solution to the normal equations for the least
squares problem
-5 -4 -3 -2 -1 0 1
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
