精通MATLAB科学计算（第3版）(王正林)10-3r

第章非线性方程求解

非线性方程是常见的一类方程，非线性方程(组)的理论远不如线性方程(组)成熟

和有效，特别是非线性方程组解的存在唯一性还没有完全解决，判断其解的存在性和解的

个数几乎没有可行的方法。利用MATLAB提供的函数，可以求解一些简单的非线性方程;

通过编程，可以解决一些较为复杂的非线性方程。

通过本章学习，读者能够熟练掌握MATLAB中的非线性方程求解相关的函数，而且

能通过编程实现多种求解非线性方程的数值算法。

10.1MATLAB中非线性方程求根函数

MATLAB中求解非线性方程的函数主要是以下两个：fzero和fsolveo其中fzero可用

来求解非线性方程，fsolve可用来求解非线性方程和非线性方程组。

10.1.1fzero函数

fzero函数用于求解单变量非线性方程1&的」基本用法及功能如下：

1.x=fzero(/；x0)求函数f在x0附近的根；

2JV=fzero(/；x0,options)使用优化选项,options可选的值很多，具体可参考MATLAB

帮助；

3.[x,fval]=fzero(/；x0,options)返回值中的fval=/(x)；

4.[x,fval,exitflag]=fzero(f,xO,options)其中exitflag标志着求解状态。

函数fzero的具体用法见下面的例子。

【例10-1]非线性方程求解函数fzero的应用实例。采用fcero函数求方程x2+x-l=0

在x=0.5附近的根。

解：在MATLAB命令窗口中输入下列命令：

»x=fzero(1xA2+x-l1,0.5)

输出计算结果为：

X=

0.6180

精通MATLAB科学计算(第2忡-----------------------------

从结果可知，方程在x=0.5附近的根为x=0.6180o

函数自er。还可以求函数在某个区间上的一个根，例如，求取方程sin(x)-0.5=0在

[0,2]区间上的根，继续在MATLAB命令窗口中输入下列命令：

»x=fzero(*sin(x)-0.5*,[02])

输出计算结果为：

X=

0.5236

从结果可知，方程在[0,2]区间上的根为x=0.5236。

需要注意的是，当用也er。求函数在某个区间上的一个根时，此函数在区间两端点处

的值符号要相反才行。

函数fzero还有一些特殊的用法，例如，求取函数f{x,y)=#-3当y=2时在x=2附

近的根，在MATLAB命令窗口中输入下列命令：

»u=inline('xAy-3','x','y');号建立非线性方程/(XJ)=*"-3

»c=fzero(uz2,[],2)

输出计算结果为：

C=

1.7321

从结果可知，所求的根为x=1.7321。

优化选项options的使用方法如下，使用0Ptimset函数，其中参数TolX选项用来设置

求解精度，在MATLAB命令窗口中输入下列命令：

»opt=optimset('TolX*,1.Oe-8);

>>x=fzero('sin(x)-0.5',0.5,opt)

x=

0.5236

>>opt=optimset(*TolX1,1.Oe-2);

>>x=fzero(*sin(x)-0.5,,0.5,opt)

x=

0.5283

从运算结果可以看出，不同精度对根有影响。

10.1.2fsolve函数

fsolve函数用于求解多变量非线性方程组。它的基本用法及功能如下：

1.x=fsolve(fun,x0)用fun定义向量函数，其定义方式可以参考例子；

2.x=fsolve(fiin,x0,options)options为优化选项；

3.[x,val]=fsolve(...)val=F(x),即函数值向量;

220►►►►

第10章非线性方程求解

4.[x,val,exitflag]=fsolve(...)exitflag标志着求解状态；

5.[x,val,exitflag,output]=fsolve(...)output包含着优化后的结果信息；

6.[x,val,exitflag,output,jacobian]=fsolve(...)jacobian为解x处的Jacobian阵。

其余参数与前面参数相似。

另外，符号数学工具箱中也有一个可用于求解方程组的函数solve，它可求各种类型方

程组的解析解。当找不到解析解时，solve会自动寻求一个近似解，并且精度很高。

需要注意的是，函数fsolve不是MATLAB符号工具箱中的函数，它位于优化工具箱

内，而solve是一个符号函数。

【例10-2】非线性方程组求解函数fsolve应用实例1。采用fsolve函数求方程组

x：_y_2=0在a,y)=(2,2)附近的根。

y2-2x-4=0

解：首先新建一个.m文件，命名为myfun.m,键入下列内容：

functionF=myfun(x)

F=[x(l)A2-x(2)-2;x(2)A2-2*x(l)-4];

在MATLAB命令窗口中输入下列命令：

»x0=[22];

>>opt=optimset(*Display','iter1);

>>[r,val]=fsolve(@myfun,xO,opt)

输出计算结果为：

NormofFirst-orderTrust-region

IterationFunc-countf(x)stepoptimalityradius

0316161

160.69900414.621

290.0001181560.1239650.03562.5

3124.3305e-0110.003045211.96e-0052.5

4156.63221e-0241.8626e-0066.78e-012

2.5

Optimizationterminated:first-orderoptimalityislessthanoptions.

TolFun.

2.21432.9032

val=

1.0e-011*

0.2075

0.1525

由计算结果可知，方程组的一个根为(x,y)=(2.2143,2.9032)。

【例10-31非线性方程组求解函数fsolve应用实例2o采用fsolve函数求矩阵方程

◄◄◄◄221

精通MATLAB科学计算（第2回

1=--的一个根，初始解向量为

24—3

解：首先新建一个.m文件，命名为myfun.m,键入下列内容：

functionF=matrixfun(x)

F=x*x-[-ll,-8;24z-3];

在MATLAB命令窗口中输入下列命令：

»xO=[l,l;lzl];

>>opt=optimset('Display',1iter');

>>r=fsolve(Gmatrixfun,xO)

输出计算结果为：

Optimizationterminated:first-orderoptimalityislessthanoptions.

TolFun.

r=

1.0000-2.0000

6.00003.0000

可以验证矩阵「一2]的平方刚好等于-11-8

24-3

10.2其他数值求根法

下面讲述的几种数值方法，在不加说明的情况下，都是在某个区间内求得方程的一个

根。其理论依据为：如果/伍）/3）<0，则在区间口力］上方程〃x）=o至少存在一个实根。

大体上可以把求解非线性方程的方法分为夹逼法（二分法和黄金分割法）和迭代法

（其他方法）两类，其中迭代法的内容十分丰富，包括一些加速和优化的技巧。下面进行

讲述。

10.2.1二分法

二分法的具体求解步骤如下。

（1）计算函数/（X）在区间⑷0中点的函数值/（等），并作蚯下面的判断：

如果/⑷/（等）V0，转到（2）；

如果>0,令,转至（1）;

如果/⑷/（岁）=0，则誓为一个根。

（2）如果a—审<£（£为预先给定的精度）,则》=券£为一个根，否则令

222►►►►

第10章非线性方程求解

6=审,转至IJ(1b

在MATLAB中编程实现的二分法的函数为：Halflntervalo

功能：用二分法求函数在某个区间上的一个零点。

调用格式：root=HalfInterval(/；a，%，eps)。

其中，/为函数名；

。为区间左端点；

6为区间右端点；

eps为根的精度；

root为求出的函数零点。

二分法的MATLAB程序代码如下：

functionroot=HalfInterval(f,a,b,eps)

考二分法求函数f在区间［a,b］上的一个零点

£函数名：f

当区间左端点；a

务区间右端点：b

软根的精度：eps

务求出的函数零点：root

if(nargin==3)

eps=l.0e-4;

end

fl=subs(sym(f),findsym(sym(f)),a);也两端点的函数值

f2=subs(sym(f),findsym(sym(f)),b);

if(fl==0)

root=a;

end

if(f2==0)

root=b;

end

if(fl*f2>0)

disp(1两端点函数值乘积大于0「)；

return;

else

root=FindRoots(f,a,b,eps);电调用求解子程序

end

functionr=FindRoots(f,a,b,eps)

f_l=subs(sym(f),findsym(sym(f)),a);

f_2=subs(sym(f),findsym(sym(f)),b);

◄◄◄◄223

精通MATLAB科学计算(第2忡-----------------------------

mf=subs(sym(f),findsym(sym(f)),(a+b)/2);%中点函数值

if(f_l*mf>0)

t=(a+b)/2;

r=FindRoots(f,t,b,eps);%右递归

else

if(f_l*mf==O)

r=(a+b)/2;

else

if(abs(b-a)<=eps)

r=(b+3*a)/4;8输出根

else

s=(a+b)/2;

r=FindRoots(f,a,s,eps);与左递归

end

end

end

【例10-4】二分法求解非线性方程应用实例。采用二分法求方程/-3x+l=0在区

间［0,1］上的一个根。

解：在MATLAB命令窗口中输入：

>>r=HalfInterval(TxA3-3*x+l',0,1)

输出计算结果为：

r=0.3473

0.3473

可见，方程d-3x+1=()在区间［0,1］上的一个根为x=0.3473。

10.2.2黄金分割法

二分法是把求解区间的长度逐次减半，而黄金分割法是把求解区间逐次缩短为前次的

0.618倍。它的求解步骤如下。

(1)设八=a+(l-0.618)*(b-a),"=0+0.618*(6-。)，且力=/("),力=/出)；

(2)如果,-司<£(给定的最小区间长度)，则输出方程的根为号■；否则转到

(3)；

(3)如果力*力<0,则令b=t2，转(1),否则如果力*/(。)>0,令a=f2,

反之令b=«,转到(1卜

在MATLAB中编程实现的黄金分割法的函数为：hjo

功能：用黄金分割法求函数在某个区间上的一个零点。

调用格式:root=hj(f,a,b,eps)。

224►►►►

第10章非线性方程求解

其中，/为函数名；

a为区间左端点；

/7为区间右端点；

eps为根的精度；

root为求出的函数零点。

黄金分割法的MATLAB程序代码如下：

functionroot=hj(f,a,b,eps)

黄金分割法求函数f在区间［a,b］上的一个零点

5t函数名：f

为区间左端点；a

8区间右端点：b

义根的精度：eps

亳求出的函数零点：root

if(nargin==3)

eps=l.Oe-4;

end

fl=subs(sym(f),findsym(sym(f)),a);

f2=subs(sym(f),findsym(sym(f)),b);

if(fl==0)

root=a;

end

if(f2==0)

root=b;

end

if(fl*f2>0)

disp(,两端点函数值乘积大于0!，)；

return;

else

tl=a+(b-a)*0.382;

t2=a+(b-a)*0.618;

f_l=subs(sym(f),findsym(sym(f)),tl);

f_2=subs(sym(f),findsym(sym(f)),t2);

tol=abs(tl-t2);

while(tol>eps)号精度控制

if(f_l*f_2<0)

a=tl;告区间两端点都调整

b=t2;

else

fa=subs(sym(f),findsym(sym(f)),a);

if(f_l*fa>0)

a=t2;扬同号左端点调整

◄◄◄◄225

精通MATLAB科学计算(第2忡--------------------------------

else

b=tl;%异号右端点调整

end

end

tl=a+(b-a)*0.382;

t2=a+(b-a)*0.618;

f_l=subs(sym(f),findsym(sym(f)),tl);

f_2=subs(sym(f),findsym(sym(f)),t2);

tol=abs(t2-tl);

end

root=(tl+t2)/2;W输出根

end

【例10-5】黄金分割法求解非线性方程应用实例。采用黄金分割法求方程

1-3x+1=0在区间［0,1］上的一个根。

解：在MATLAB命令窗口中输入：

»r=hj('xA3-3*x+l',0,l)

输出计算结果为：

r=0.3473

可见，方程在区间［0,1］上的一个根为x=0.3473,与二分法求出来的结果相同。

10.2.3不动点迭代法

求方程/(x)=0的根，可以先把方程改写成如下的形式：

x=/(x)+x

于是得到不动点迭代法的其中一种迭代公式：

X„=/(x„-i)+x„_1

此种不动点迭代法很有可能不收敛，因为它的本质是求函数y=/(x)+x与直线y=x

的交点，而它们不一定存在交点。即使收敛，其速度也十分慢，因此有了艾特肯加速迭代

与史蒂芬森加速迭代。

在MATLAB中编程实现的不动点迭代法的函数为：StablePointo

功能：用不动点迭代法求函数的一个零点。

调用格式：［root,n］=StablePoint(f,x0,eps)。

其中，/为函数名；

M)为初始迭代向量；

eps为根的精度；

226►►►►

第10章非线性方程求解

root为求出的函数零点；

〃为迭代步数。

不动点迭代法的MATLAB程序代码如下：

function[root,n]=StablePoint(f,xO,eps)

常用不动点迭代法求函数的一个零点

3初始迭代向量：x0

卷根的精度：eps

%求出的函数零点：root

%迭代步数：n

if(nargin==2)

eps=l.Oe-4;

end

tol=l;

root=x0;

n=0;

while(tol>eps)

n=n+l;

rl=root;

root=subs(sym(f),findsym(sym(f)),rl)+rl;超迭代的核心公式

tol=abs(root-rl);

end

【例10-6】不动点迭代法求解非线性方程应用实例。采用不动点迭代法求方程

-y=+x-2=0的一个根。

解：在MATLAB命令窗口中输入：

>>[rzn]=StablePoint(11/sqrt(x)+x-2*,0.5)

r=0.3820

nV.J7QJJon2

n=4

二

从计算结果可以看出，经过4步迭代，得出方程J=+x-2=0的一个根为

x=0.3820o

1.艾特肯加速迭代法

艾特肯加速迭代法是在计算出X〃,X〃+l,X〃+2后，对作以下修正：

*(X〃+1—x”)~

x"+i=-

x“+2—2x〃+i+xn

Y◄◄◄227

精通MATLAB科学计算(第2忡-----------------------------

然后用焉+i来逼近方程的根。

艾特肯加速迭代法是先用不动点迭代法算出一系列{4},再对此系列作修正得到系列

{xn},用后者来逼近方程的根。

在MATLAB中编程实现的艾特肯加速迭代法的函数为：AtkenStablePointo

功能：用艾特肯加速迭代法求函数的一个零点。

调用格式:[root,/?]=AtkenStablePoint(f,xO,eps)o

其中，/为函数名；

M)为初始迭代向量；

eps为根的精度；

root为求出的函数零点；

〃为迭代步数。

艾特肯加速迭代法的MATLAB程序代码如下：

function[root,n]=AtkenStablePoint(f,xO,eps)

Z用艾特肯加速迭代法求函数f的一个零点

常初始迭代向量：xO

Z根的精度：eps

年求出的函数零点：r。”

金迭代步数：n

if(nargin==2)

eps=l.Oe-4;

end

tol=l;

root=xO;

x(1:2)=0;

n=0;

m=0;

a2=x0;

while(tol>eps)

n=n+l;

al=a2;

rl=root;

root=subs(sym(f),findsym(sym(f)),rl)+rl;%求出根

x(n)=root;w把所有的根都保存下来

if(n>2)

m=m+1;

a2=x(m)-(x(m+1)-x(m))^2/(x(m+2)-2*x(m+1)+x(m));%对根进行艾特肯修正

tol=abs(a2-al);

end

end

root=a2;

228►►►►

-------------------------------------------------------第10章非线性方程求解

【例10-7】艾特肯加速迭代法求解非线性方程应用实例。采用艾特肯加速迭代法求

方程的一个根，迭代初始值为。

耳+x-2=00.5

解：在MATLAB命令窗口中输入：

>>[r,n]=AtkenStablePoint(11/sqrt(x)+x-2,,0.5)

r=0.3820

n=4

-4

从计算结果可以看出，经过4步迭代，得出方程9+x-2=0的一个根为x=0.3820。

这个方程比较简单，所以通过简单的几步迭代就可得出根。从下面的结果可以看出艾特肯

加速迭代法的优点。

»[rzn]=AtkenStablePoint(*1/sqrt(x)+x-210.999)

r=1.0000

Ji-•n7o2n72A

n=£

4

>>[r,n]=StablePoint('1/sqrt(x)+x-210.999)

r=0.3820

0.3820

n=21

w

上面的例子说明采用初始值0.999时，经过4步艾特肯加速迭代法，可以得到方程

+x-2=0的另一个根x=l,而普通的不动点迭代法却得不到。

Y◄◄◄229

精通MATLAB科学计算(第2忡-----------------------------

2.史蒂芬森加速迭代法

史蒂芬森加速与艾特肯加速不同的地方在于前者是在迭代的同时就进行修正，它的迭

代公式为：

yn=/(x“)+X"

z”=/仇)+%

(%-X")?

X„+l=X„------------------

z“-2Ml+xn

在MATLAB中编程实现的史蒂芬森加速迭代法的函数为：StevenStablePointo

功能：用史蒂芬森加速迭代法求函数的一个零点。

调用格式：[root,«]=StevenStablePoint(/；.rO,eps)o

其中，/为函数名；

xO为初始迭代向量；

eps为根的精度；

root为求出的函数零点；

〃为迭代步数。

史蒂芬森加速的MATLAB程序代码如下：

function[root,n]=StevenStablePoint(f,xO,eps)

徒用史蒂芬森加速迭代法求函数f的一个零点

%初始迭代向量：xO

%根的精度：eps

学求出的函数零点：root

为迭代步数：n

if(nargin==2)

eps=l.Oe-4;

end

tol=l;

root=xO;

n=0;

while(tol>eps)

n=n+l;

rl=root;

y=subs(sym(f),findsym(sym(f)),rl)+rl;

z=subs(sym(f),findsym(sym(f)),y)+y;

root=rl-(y-rl)A2/(z-2*y+rl);%对每次算出的根立即进行修正

tol=abs(root-rl);

end

【例10-8】史蒂芬森加速迭代法求解非线性方程应用实例。采用史蒂芬森加速迭代

230►►►►

第10章非线性方程求解

法求方程9+x-2=0的一个根。

解：在MATLAB命令窗口中输入：

>>[r,n]=StevenStablePoint('1/sqrt(x)+x-2',0.999)

r=1.0000

一i•Vn>nVoVoV

n=2

2

>>[r,n]=StevenStablePoint(11/sqrt(x)+x-2'z0.5)

r=0.3820

n=4

4

1

对比上例可以看出，史蒂芬森加速迭代法不仅能求出方程

忑+x-2=0的一个根

x=0.3820,而且它比艾特肯加速迭代法更快地求出另一个根x=l。

10.2.4弦截法

弦截法的算法过程如下：

(1)过两点(a,/(a)),(6J3))作一直线，它与X轴有一个交点，记为汨；

(2)如果f[a}f{x\)<0,过两点(aj(a)),(x1,/(的))作一直线，它与X轴的交点记

为,否则过两点(6,7(b)),(X[,/但))作一直线，它与x轴的交点记为打；

(3)如此下去，直到,就可以认为X“为/(尤)=0在区间[a,b]上的一•/

根。

xi-a

Xk=a-/(«)/W(xi)<0

(4)为,的递推公式为：.

Xk-i-b

Xk=b-/3)/(a)/(x_,)>0

/(%*._I)-/(6)t

且Xi=a----------------/(a)o

在MATLAB中编程实现的弦截法的函数为：Secanto

功能：用弦截法求函数在某个区间上的一个零点。

调用格式:root=Secant(f,a,h,eps)o

其中，/为函数名；

a为区间左端点；

6为区间右端点；

eps为根的精度；

Y◄◄◄231

精通MATLAB科学计算(第2贿——

root为求出的函数零点。

弦截法的MATLAB程序代码如下：

functionroot=Secant(fza,b,eps)

%弦截法求函数f在区间上的一个零点

£函数名：f

他区间左端点；a

%区间右端点：b

为根的精度：eps

星求出的函数零点：root

if(nargin==3)

eps=l.Oe-4;

end

fl=subs(sym(f),findsym(sym(f)),a);

f2=subs(sym(f),findsym(sym(f)),b);

if(fl==0)

root=a;

end

if(f2==0)

root=b;

end

if(fl*f2>0)

disp('两端点函数值乘积大于01;

return;

else

tol=l;

fa=subs(sym(f),findsym(sym(f)),a);

fb=subs(sym(f),findsym(sym(f)),b);

root=a-(b-a)*fa/(fb-fa);多迭代初始值

while(tol>eps)

rl=root;

fx=subs(sym(f),findsym(sym(f)),rl);

s=fx*fa;

if(s==0)

root=rl;

else

if(s>0)

root=b-(rl-b)*fb/(fx-fb);

else

root=a-(rl-a)*fa/(fx-fa);

end

end

tol=abs(root-rl);

end

end

232►►►>

-------------------------------------------------------第10章非线性方程求解

【例10-9】弦截法求解非线性方程应用实例。采用弦截法求方程lgx+&=2在区

间［1,4］上的一个根。

解：在MATLAB命令窗口中输入：

1

>>r=Secant(*log(x)+sqrt(x)-2,lz4)

输出计算结果为：

r=1.8773

1.8773

由计算结果可知，lgx+&=2在区间［1,4］上的一个根为x=1.8773。

10.2.5史蒂芬森弦截法

史蒂芬森弦截法是弦截法的一种变形，它的递推公式为：

Xk=Xk-\---------------'"-------------f(XA-1)

/(Xz+fg-)八

且有

x\=a--------------------------f(a)o

史蒂芬森弦截法比弦截法的精度要高，收敛速度要快。

在MATLAB中编程实现的史蒂芬森弦截法的函数为：StevenSecanto

功能：用史蒂芬森弦截法求函数在某个区间上的一个零点。

调用格式:root=StevenSecant(f,a,b,eps)o

其中，/为函数名；

a为区间左端点；

6为区间右端点；

eps为根的精度；

root为求出的函数零点。

史蒂芬森弦截法的MATLAB程序代码如下：

functionroot=StevenSecant(f,a,b,eps)

亳史蒂芬森弦截法求函数f在区间［a,b］上的一个零点

会函数名：f

》区间左端点；a

彩区间右端点：b

国根的精度：eps

省求出的函数零点：root

Y◄◄◄233

精通MATLAB科学计算(第2忡-----------------------------

if(nargin==3)

eps=l.Oe-4;

end

fl=subs(sym(f),findsym(sym(f))za);

f2=subs(sym(f),findsym(sym(f)),b);

if(fl==O)

root=a;

end

if(f2==0)

root=b;

end

if(fl*f2>0)

dispC两端点函数值乘积大于0!，)；

return;

else

tol=l;

fa=subs(sym(f),findsym(sym(f)),a);

fb=subs(sym(f),findsym(sym(f)),b);

faa=subs(sym(f),findsym(sym(f)),fa+a);

root=a-fa*fa/(faa-fa);多迭代初始值

while(tol>eps)

rl=root;

fx=subs(sym(f),findsym(sym(f)),rl);

v=fx+rl;

fxx=subs(sym(f),findsym(sym(f)),v);

root=rl-fx*fx/(fxx-fx);S递推的核心公式

tol=abs(root-rl);

end

end

【例10-10］史蒂芬森弦截法求解非线性方程应用实例。采用史蒂芬森弦截法求方

程lgx+=2在区间［1,4］上的一个根。

解：在MATLAB命令窗口中输入：

>>r=StevenSecant(*log(x)+sqrt(x)-2',1,4)

输出计算结果为：

r=1

一工

结果出现错误，改正方法可以是把求解区间缩短，缩短为［1.5,2］：

>>r=StevenSecant(*log(x)+sqrt(x)-2',1.5,2)

r=1.8773

-3

这样就得出了正确的根1.8773。从上面的例子可以看出求解区间的选择也是很重要的，

一般来说，求解区间越短越好。

234►►►►

-------------------------------------------------------第；0章非线性方程求解

10.2.6抛物线法

弦截法其实是用不断缩短的直线来近似函数/（X）,而抛物线法则采用三个点来近似函

数/（X），并且用抛物线与横轴的交点来逼近函数./a）的根。它的算法过程如下：

（1）选定初始值工0,两,切,并计算/（%0）,/（两）,/（x2）和以下差分:

即-Xo

/1电的,知=/8闻一/氏项

X2-XQ

一般取XO=Q,x\=b,a<X2<bo注意不要使三点共线。

（2）用牛顿插值法对三点（工0,/'（工0））,（两,/（两））,（工2J（X2））进行插值得到一条抛物线,

它有两个根：

-g±Vg2-4AC

%3=+

2C

其中

A=f(x2),C=f[X2,X\,Xo],

B=f[x2,X1]+f[x2,X|,Xo](x2-XI)。

两个根中只取靠近X2的那个根，即土号取与8同号，

即

2A

X3=X2------------------ff

B+sgn(B)ylB2-4AC

（3）用百户2,看代替X。,即,X2,重复以上步骤，并有以下递推公式：

2An

B,i+sgn(^z/)\JB,r-4AnCn

其中

A”=,/(x"),C”=f[xn,Xn-\,Xn-2],

Bn=f[xn,Xn-\]+f[xnXn-2](x„-Xn-\)»

（4）进行精度控制。

在MATLAB中编程实现的抛物线法的函数为：Parabolao

功能：用抛物线法求函数在某个区间上的一个零点。

调用格式：root=Parabola（f,a,h,x,eps）o

词<<<235

精通MATLAB科学计算(第2忡-----

其中，/为函数名；

〃为区间左端点；

人为区间右端点；

X为初始迭代点；

eps为根的精度；

root为求出的函数零点。

抛物线法的MATLAB程序代码如下：

functionroot=Parabola(f,a,b,x,eps)

乡抛物线法求函数f在区间［a,b］上的一个零点

W函数名：f

阳区间左端点；a

务区间右端点：b

》初始迭代点：x

他根的精度：eps

%求出的函数零点：root

if(nargin==4)

eps=l.Oe-4;

end

fl=subs(sym(f),findsym(sym(f)),a);

f2=subs(sym(f),findsym(sym(f)),b);

if(fl==0)

root=a;

end

if(f2==0)

root=b;

end

if(fl*f2>0)

disp('两端点函数值乘积大于0!，)；

return;

else

tol=l;

fa=subs(sym(f),findsym(sym(f)),a);

fb=subs(sym(f),findsym(sym(f)),b);

fx=subs(sym(f),findsym(sym(f)),x);

dl=(fb-fa)/(b-a);

d2=(fx-fb)/(x-b);

d3=(d2-dl)/(x-a);

B=d2+d3*(x-b);

root=x-2*fx/(B+sign(B)*sqrt(BA2-4*fx*d3)

t=zeros(3);

t(1)=a;

236►►►►

第10章非线性方程求解

t(2)=b;

t(3)=x;

while(tol>eps)

t(l)=t(2);称保存三个点

t(2)=t(3);

t(3)=root;

fl=subs(sym(f),findsym(sym(f)),t(1));%计算三点的函数值

f2=subs(sym(f),findsym(sym(f))zt(2));

f3=subs(sym(f),findsym(sym(f)),t⑶)；

dl=(f2-fl)/(t(2)-t(l));告计算三个差分

d2=(f3-f2)/(t(3)-t(2));

d3=(d2-dl)/(t(3)-t(l));

B=d2+d3*(t(3)-t(2));%计算算法中的B

root=t(3)-2*f3/(B+sign(B)*sqrt(BA2-4*f3*d3));

tol=abs(root-t(3));

end

end

【例10-11】抛物线法求解非线性方程应用实例I。采用抛物线法求方程

lgx+4=2在区间［1,4］上的一个根，迭代初始值为2。

解：在MATLAB命令窗口中输入：

>>“Parabola('sqrt(x)+log(x)-2',1,4,2)

r=1.8773

1.8773

【例10-12】抛物线法求解非线性方程应用实例2o分别从两个区间［0.5,1.5］和［0.1,1］

1

上采用抛物线法求方程的根。

忑+x-2=0

解：在MATLAB命令窗口中输入：

»r=Parabola('1/sqrt(x)+x-2\0.5,1.5,0.8)告迭代区间为[0・5,1.5],初始值为0.8

r=1

工

>>r=Parabola(*1/sqrt(x)+x-2*z0.1zl,0.5)带迭代区间为[0.1,1],初始值为0.5

r=0.3820

1

由此可以看出抛物线法也能求出方程

忑+x-2=0的两个根。

10.2.7牛顿法

弦截法本质上是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一

种切线法，它从一端向一个方向逼近方程的根，其递推公式为：

Y◄◄◄237

精通MATLAB科学计算(第2回

fM

初始值可以取/(幻和/(b)的较大者，这样可以加快收敛速度。

和牛顿法有关的还有简化牛顿法和牛顿下山法。

在MATLAB中编程实现的牛顿法的函数为：NewtonRooto

功能：用牛顿法求函数在某个区间上的一个零点。

调用格式：root=NewtonRoot(f,a,b,eps)。

其中，/为函数名；

。为区间左端点；

6为区间右端点；

eps为根的精度；

root为求出的函数零点。

牛顿法的MATLAB程序代码如下：

functionroot=NewtonRoot(f,a,b,eps)

W牛顿法求函数f在区间［a,b］上的一^零点

专函数名：f

得区间左端点；a

专区间右端点：b

传根的精度：eps

务求出的函数零点：root

if(nargin==3)

eps=l.Oe-4;

end

fl=subs(sym(f),findsym(sym(f)),a);

f2=subs(sym(f),findsym(sym(f)),b);

if(fl==0)

root=a;

end

if(f2==0)

root=b;

end

if(fl*f2>0)

disp(1两端点函数值乘积大于0!，)；

return;

else

tol=l;

fun=diff(sym(f));%求导数

fa=subs(sym(f),findsym(sym(f)),a);

238►►►►

第10章非线性方程求解

fb=subs(sym(f),findsym(sym(f)),b);

dfa=subs(sym(fun),findsym(sym(fun)),a);

dfb=subs(sym(fun),findsym(sym(fun)),b);

if(dfa>dfb)考初始值取两端点导数较大者

root=a-fa/dfa;

else

root=b-fb/dfb;

end

while(tol>eps)

rl=root;

fx=subs(sym(f),findsym(sym(f)),rl);

dfx=subs(sym(fun),findsym(sym(fun)),rl);为求该点的导数值

root=rl-fx/dfx;告迭代