Matlab Math
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
零点和根
Bisection
我们知道 在 1和 2之间，先试 x=3/2,因为 x的平
方大于 2，那么 在 [1,3/2]…
2
2
,...64271,32131,1651,831,411,211
M=2;a=1;b=3;k=0;
while b-a>eps
x=(a+b) /2;
if x^2 >M
b=x
else
a = x
end
k=k+1;
end
Interval bisection
观察一下 a,b变化的情况：
程序结束需要多少次？
最后 a,b是多少？
用 format hex表示呢？
为什么是这样的？
Interval bisection－－ f(x)
Find a very small interval,perhaps two successive floating
point numbers,on which the function change sign.
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;
end
f = inline('x^2-2'); a=1; b=2; 计算 sqrt(2)
Interval bisection
二分法是慢的，用上面的代码，对任何
函数她都需要运行 52步，但是她是完全
可靠的。这与我们后面介绍的方法不同。
Newton's Method
对于求 f(x) = 0的根,Newton法要求过 y=f(x)上任一点
的切线，并求切线与 x-轴的交。迭代过程是
)('
)(
1
n
n
nn xf
xfxx ??
?
k=0;
while abs(x-xprev)>eps*abs(x)
xprev = x;
x = x-f(x)/fprime(x)
end
Newton's Method
对于平方根问题,Newton法是优雅而有效的
Mxxf ?? 2)( M零点就是
???
?
???
? ?????
?
n
n
n
n
nn x
Mx
x
Mxxx
2
1
2
2
1
Newton法
k=0;
while abs(x-xprev)>eps*abs(x)
xprev = x;
x = 0.5*(x+M/x)
end
1.50000000000000
1.41666666666667
1.41421568627451
1.41421356237469
1.41421356237309
1.41421356237309
6步
Newton's Method
缺点：
函数 f(x)必须是光滑的，要容易计算导数 ;初始猜想要靠近
零点。
优点：收敛速度快，二次收敛。如果函数的一阶和二阶导数
连续，且初始值靠近零点，我们可以证明
2
1 )('
)(''
2
1
n
n
n exf
fe ???
? )( 21 nn eOe ??
Newton's Method
)('
)(
1
n
n
nn xf
xfxx ??
?
下面构造一个让 Newton方法无限进行而不收敛的例子
)(1 axax nn ????? )(
)('
)( ax
xf
xfax ?????
)(2
1
)(
)('
axxf
xf
?? ||)()( axaxs i g nxf ???
Newton方法不收敛，原因？
0 0, 5 1 1, 5 2 2, 5 3 3, 5 4
- 1, 5
-1
- 0, 5
0
0, 5
1
1, 5
x
s i g n ( x - 2 ) s q r t ( abs ( x - 2 ))
Secant Method
用有限差分来代替导数的计算
1
1 )()(
?
?
?
??
nn
nn
n xx
xfxfs
n
n
nn s
xfxx )(
1 ???
)('
)(
1
n
n
nn xf
xfxx ??
?
k=0;
while abs(b-a)>eps*abs(b)
c=a;a=b;
b=b+(b-c)/(f(c)/f(b)-1);
k=k+1;
end
1.33333333333333
1.40000000000000
1.41463414634146
1.41421143847487
1.41421356205732
1.41421356237310
1.41421356237310
7步
13
1
1 )('
)(')(')(''
2
1
?
?
? ?? nn
nn
n eef
fffe
?
???
Secant Method
如果 f(x)的一阶导数和二阶导数连续
)( 11 ?? ? nnn eeOe Superline convergence
2
51),(
1
???
? ?
?
nn eOe
即每作一次迭代有效位将增加 1.6位
再考虑一下 ||)()( axaxs i g nxf ???
Inverse Quadratic Interpolation
Secant方法用前两个点来确定第三个点，能不能用三个点呢？
))(() ),(() ),(( cfPcbfPbafPa ??? IQI
k=0;
while abs(c-b)>eps*abs(c)
x=polyinterp([f(a),f(b),f(c)],[a,b,c],0);
a=b;b=c;c=x;
k=k+1;
end
1.66666666666667
1.40151515151515
1.41389545042106
1.41421378833549
1.41421356237333
1.41421356237310
1.41421356237309
[a,b,c]=0,1,2
Inverse Quadratic Interpolation
a,b,c = -2,0,2
a,b,c = -2.00001,0,1.99999
看一下下面两组数据,IQI方法会发生什么情况？
Zeroin
? Start with a and b so that f(a) and f(b) have opposite signs
? Use a secant step to give c between a and b
?Repeat the following steps until |b-a|<eps|b| or f(b)=0
?Arrange a,b and c so that
? f(a) and f(b) have opposite signs
? |f(b)|<=|f(a)|
? c is the previous value of b
? If c~=a,consider an IQI step
?If c=a,consider a secant stop
?If the IQI or secant step is in the interval [a,b],take it
? If the step is not in the interval,use bisection.
fzerotx
Matlab中 zeroin的算法实现是 fzero.
x = fzero(FUN,x0) %x0可以是数，或区间
x = fzero(‘besselj(0,x)’,0) %?
x = fzero(‘besselj(0,x)’,[0,pi])
x = fzerotx(‘besselj(0,x)’,[0,pi])
- 3 0 - 2 0 - 1 0 0 10 20 30
- 0, 5
0
0, 5
1
x
b e s s e l j ( 0,x )
ezplot(‘besselj(0,x)’,[-30,30])
)(0 xJ
函数调用
?字符串表达式 ‘ cos(pi*t)’ ‘besselj(0,x)’
? Inline函数
?Function handle @bessj0
z = fzerotx(@bessj0,[0,pi])
function y=bessj0(x)
y=besselj(0,x)
?符号表示 syms x
F=besselj(0,x)
z=fzerotx(F,[0,pi])
F=inline(‘besselj(0,x)-y’,’x’,’y’)
x = fzerotx(F,[0,2],.5)
fzerogui
zeroin 的图形化描述
fzerogui(‘besselj(0,x)’,[0 3.83]);
Fzerogui.m
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
零点和根
Bisection
我们知道 在 1和 2之间，先试 x=3/2,因为 x的平
方大于 2，那么 在 [1,3/2]…
2
2
,...64271,32131,1651,831,411,211
M=2;a=1;b=3;k=0;
while b-a>eps
x=(a+b) /2;
if x^2 >M
b=x
else
a = x
end
k=k+1;
end
Interval bisection
观察一下 a,b变化的情况：
程序结束需要多少次？
最后 a,b是多少？
用 format hex表示呢？
为什么是这样的？
Interval bisection－－ f(x)
Find a very small interval,perhaps two successive floating
point numbers,on which the function change sign.
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;
end
f = inline('x^2-2'); a=1; b=2; 计算 sqrt(2)
Interval bisection
二分法是慢的，用上面的代码，对任何
函数她都需要运行 52步，但是她是完全
可靠的。这与我们后面介绍的方法不同。
Newton's Method
对于求 f(x) = 0的根,Newton法要求过 y=f(x)上任一点
的切线，并求切线与 x-轴的交。迭代过程是
)('
)(
1
n
n
nn xf
xfxx ??
?
k=0;
while abs(x-xprev)>eps*abs(x)
xprev = x;
x = x-f(x)/fprime(x)
end
Newton's Method
对于平方根问题,Newton法是优雅而有效的
Mxxf ?? 2)( M零点就是
???
?
???
? ?????
?
n
n
n
n
nn x
Mx
x
Mxxx
2
1
2
2
1
Newton法
k=0;
while abs(x-xprev)>eps*abs(x)
xprev = x;
x = 0.5*(x+M/x)
end
1.50000000000000
1.41666666666667
1.41421568627451
1.41421356237469
1.41421356237309
1.41421356237309
6步
Newton's Method
缺点：
函数 f(x)必须是光滑的，要容易计算导数 ;初始猜想要靠近
零点。
优点：收敛速度快，二次收敛。如果函数的一阶和二阶导数
连续，且初始值靠近零点，我们可以证明
2
1 )('
)(''
2
1
n
n
n exf
fe ???
? )( 21 nn eOe ??
Newton's Method
)('
)(
1
n
n
nn xf
xfxx ??
?
下面构造一个让 Newton方法无限进行而不收敛的例子
)(1 axax nn ????? )(
)('
)( ax
xf
xfax ?????
)(2
1
)(
)('
axxf
xf
?? ||)()( axaxs i g nxf ???
Newton方法不收敛，原因？
0 0, 5 1 1, 5 2 2, 5 3 3, 5 4
- 1, 5
-1
- 0, 5
0
0, 5
1
1, 5
x
s i g n ( x - 2 ) s q r t ( abs ( x - 2 ))
Secant Method
用有限差分来代替导数的计算
1
1 )()(
?
?
?
??
nn
nn
n xx
xfxfs
n
n
nn s
xfxx )(
1 ???
)('
)(
1
n
n
nn xf
xfxx ??
?
k=0;
while abs(b-a)>eps*abs(b)
c=a;a=b;
b=b+(b-c)/(f(c)/f(b)-1);
k=k+1;
end
1.33333333333333
1.40000000000000
1.41463414634146
1.41421143847487
1.41421356205732
1.41421356237310
1.41421356237310
7步
13
1
1 )('
)(')(')(''
2
1
?
?
? ?? nn
nn
n eef
fffe
?
???
Secant Method
如果 f(x)的一阶导数和二阶导数连续
)( 11 ?? ? nnn eeOe Superline convergence
2
51),(
1
???
? ?
?
nn eOe
即每作一次迭代有效位将增加 1.6位
再考虑一下 ||)()( axaxs i g nxf ???
Inverse Quadratic Interpolation
Secant方法用前两个点来确定第三个点，能不能用三个点呢？
))(() ),(() ),(( cfPcbfPbafPa ??? IQI
k=0;
while abs(c-b)>eps*abs(c)
x=polyinterp([f(a),f(b),f(c)],[a,b,c],0);
a=b;b=c;c=x;
k=k+1;
end
1.66666666666667
1.40151515151515
1.41389545042106
1.41421378833549
1.41421356237333
1.41421356237310
1.41421356237309
[a,b,c]=0,1,2
Inverse Quadratic Interpolation
a,b,c = -2,0,2
a,b,c = -2.00001,0,1.99999
看一下下面两组数据,IQI方法会发生什么情况？
Zeroin
? Start with a and b so that f(a) and f(b) have opposite signs
? Use a secant step to give c between a and b
?Repeat the following steps until |b-a|<eps|b| or f(b)=0
?Arrange a,b and c so that
? f(a) and f(b) have opposite signs
? |f(b)|<=|f(a)|
? c is the previous value of b
? If c~=a,consider an IQI step
?If c=a,consider a secant stop
?If the IQI or secant step is in the interval [a,b],take it
? If the step is not in the interval,use bisection.
fzerotx
Matlab中 zeroin的算法实现是 fzero.
x = fzero(FUN,x0) %x0可以是数，或区间
x = fzero(‘besselj(0,x)’,0) %?
x = fzero(‘besselj(0,x)’,[0,pi])
x = fzerotx(‘besselj(0,x)’,[0,pi])
- 3 0 - 2 0 - 1 0 0 10 20 30
- 0, 5
0
0, 5
1
x
b e s s e l j ( 0,x )
ezplot(‘besselj(0,x)’,[-30,30])
)(0 xJ
函数调用
?字符串表达式 ‘ cos(pi*t)’ ‘besselj(0,x)’
? Inline函数
?Function handle @bessj0
z = fzerotx(@bessj0,[0,pi])
function y=bessj0(x)
y=besselj(0,x)
?符号表示 syms x
F=besselj(0,x)
z=fzerotx(F,[0,pi])
F=inline(‘besselj(0,x)-y’,’x’,’y’)
x = fzerotx(F,[0,2],.5)
fzerogui
zeroin 的图形化描述
fzerogui(‘besselj(0,x)’,[0 3.83]);
Fzerogui.m
