第二章非线性方程(组)的数值解法的matlab程序

本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（M üller ）法和迭代法的加速等）及其MATLAB 程序，求解非线性方程组的方法及其MATLAB 程序.

2.1 方程（组）的根及其MATLAB 命令

2.1.2 求解方程（组）的solve 命令

求方程f (x )=q (x )的根可以用MATLAB 命令:

>> x=solve('方程f(x)=q(x)',’待求符号变量x ’)

求方程组f i (x 1,…,x n )=q i (x 1,…,x n ) (i =1,2,…,n )的根可以用MATLAB 命令:

>>E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)'); ……………………………………………………. En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)'); [x1,x2,…,xn]=solve(E1,E2,…,En, x1,…,xn)

2.1.3 求解多项式方程（组）的roots 命令

如果)(x f 为多项式，则可分别用如下命令求方程0)(=x f 的根，或求导数)('x f （见表 2-1）.

2.1.4 求解方程（组）的fsolve 命令

如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots 命令.如果非线性方程（组）是含有超越函数，则无法使用roots 命令，需要调用MATLAB 系统中提供的另一个程序fsolve 来求解.当然，程序fsolve 也可以用于多项式方程（组），但是它的计算量明显比roots 命令的大.

fsolve 命令使用最小二乘法（least squares method ）解非线性方程（组）

(F X =)0

的数值解，其中X 和F （X ）可以是向量或矩阵.此种方法需要尝试着输入解X 的初始值（向量或矩阵）X 0，即使程序中的迭代序列收敛，也不一定收敛到(F X =)0的根（见例2.1.8）.

fsolve 的调用格式： X=fsolve(F,X0)

输入函数)(x F 的M 文件名和解X 的初始值（向量或矩阵）X 0，尝试着解方程（组）

(F X =)0，运行后输出(F X =)0解的估计值（向量或矩阵）X .

要了解更多的调用格式和功能请输入:help fsolve ,查看说明.

2.2 搜索根的方法及其MATLAB 程序

求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根，需要搜索根所在的区间或确定根的初始近似值（简称初始值）.搜索根的近似位置的常用方法有三种:作图法、逐步搜索法和二分法等，使用这些方法的前提是高等数学中的零点定理.

2.2.1 作图法及其MATLAB 程序

作函数的图形的方法很多，如用计算机软件的图形功能画图，或用高等数学中应用导数作图，或用初等数学的函数叠加法作图等.下面介绍两种作图程序.

作函数)(x f y =在区间 [a,b ] 的图形的MATLAB 程序一

x=a:h:b; % h 是步长 y=f(x); plot(x,y) grid, gtext('y=f(x)')

说明：⑴ 此程序在MATLAB 的工作区输入，运行后即可出现函数)(x f y =的图形.此图形与x 轴交点的横坐标即为所要求的根的近似值.

⑵ 区间[a,b ] 的两个端点的距离 b-a 和步长h 的绝对值越小，图形越精确. 作函数)(x f y =在区间 [a,b ]上的图形的MATLAB 程序二 将)(x f y =化为)()(x g x h =，其中)()(x g x h 和是两个相等的简单函数

x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2)

grid,gtext(' y1=h(x),y2=g(x)') 说明：此程序在MATLAB 的工作区输入，运行后即可出现函数)()(21x g y x h y ==和的图形.两图形交点的横坐标即为所要求的根的近似值.

下面举例说明如何用计算机软件MATLAB 的图形功能作图.

2.2.2 逐步搜索法及其MATLAB 程序

逐步搜索法也称试算法.它是求方程0)(=x f 根的近似值位置的一种常用的方法. 逐步搜索法依赖于寻找连续函数)(x f 满足)(a f 与)(b f 异号的区间],[b a .一旦找到区间，无论区间多大，通过某种方法总会找到一个根.

MATLAB 的库函数中没有逐步搜索法的程序，现提供根据逐步搜索法的计算步骤和它的

% 输入的量--- a 和b 是闭区间[a,b]的左、右端点； %---h 是步长；

%---tol 是预先给定的精度.

% 运行后输出的量---k 是搜索点的个数；

% --- r 是方程 在[a,b]上的实根的近似值，其精度是tol ； X=a:h:b;Y=funs(X);n=(b-a)/h+1;m=0; X(n+1)=X(n);Y(n+1)=Y(n); for k=2:n

X(k)=a+k*h;Y(k)=funs(X(k)); %程序中调用的funs.m 为函数

sk=Y(k)*Y(k-1); if sk<=0,

m=m+1;r(m)=X(k); end

xielv=(Y(k+1)-Y(k))*(Y(k)-Y(k-1)); if (abs(Y(k))

例2.2.4 用逐步搜索法的MATLAB 程序分别求方程0332223=--+x x x 和

0)2sin(cos 3=x 在区间]2,2[-上的根的近似值，要求精度是0.000 1.

解 在MATLAB 工作窗口输入如下程序

>> [k,r]=zhubuss(-2,2,0.001,0.0001)

运行后输出的结果

k =4001

r = -1.2240 -1.0000 -1.0000 -0.9990 1.2250

即搜索点的个数为k =4 001，其中有5个是方程⑴的近似根，即r = -1.224 0，-1.000 0，-1.000 0，-0.999 0，1.225 0，其精度为0.000 1.

在程序中将y =2.*x .^3+2.*x .^2-3.*x -3用y =sin(cos(2.*x .^3)) 代替，可得到方程⑵在区间]2,2[-上的根的近似值如下

r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.9230 1.3310

1.5780 1.7650 1.9200

如果读者分别将方程⑴的结果与例2.2.3比较，方程⑵的结果与例2.1.2比较，将会发现逐步搜索法的MATLAB 程序的优点.如果精度要求比较高，用这种逐步搜索法是不合算的.

2.3 二分法及其MATLAB 程序

2.3.1 二分法

二分法也称逐次分半法.它的基本思想是：先确定方程0)(=x f 含根的区间（a,b ），再把区间逐次二等分.

我们可以根据式（2.3b ）、区间[a,b ]和误差ε，编写二分法求方程根的迭代次数的MATLAB

2.3.2 二分法的MATLAB 程序

二分法需自行编制程序，现提供用二分法求方程f (x )=0的根*

x 的近似值k x 的步骤和式

a(1)=a; b(1)=b;

ya=fun(a(1)); yb=fun(b(1)); %程序中调用的fun.m 为函数 if ya* yb>0,

disp('注意：ya*yb>0,请重新调整区间端点a 和b.'), return end

max1=-1+ceil((log(b-a)- log(abtol))/ log(2)); % ceil 是向∞+ 方向取整 for k=1: max1+1

a;ya=fun(a); b;yb=fun(b); x=(a+b)/2; yx=fun(x); wuca=abs(b-a)/2; k=k-1; [k,a,b,x,wuca,ya,yb,yx] if yx==0

a=x; b=x; elseif yb*yx>0

b=x;yb=yx; else

a=x; ya=yx; end

if b-a< abtol , return , end end

k=max1; x; wuca; yx=fun(x);

例2.3.3 确定方程x 3-x +4=0的实根的分布情况，并用二分法求在开区间 (-2,-1)内的实根的近似值，要求精度为0.001.

解 （1）先用两种方法确定方程x 3-x +4=0 的实根的分布情况。

方法1 作图法.

在MATLAB 工作窗口输入如下程序

>>x=-4:0.1:4;

y=x.^3-x +4; plot(x,y) grid,gtext('y=x^3-x+4')

画出函数f (x )=x 3-x+4的图像，如图2-7.从图像可以看出，此曲线有两个驻点3

3

±

都在x 轴的上方，在(-2,-1)内曲线与x 轴只有一个交点，则该方程有唯一一个实根，且在 (-2,-1)内.

方法2 试算法.

在MATLAB 工作窗口输入程序

>>x=-4: 1:4,y=x.^3-x +4

运行后输出结果

x = -4 -3 -2 -1 0 1 2 3 4

y = -56 -20 -2 4 4 4 10 28 64

由于连续函数f (x )满足0)1()2(<-?-f f ，所以此方程在(-2,-1)内有一个实根.

（2） 用二分法的主程序计算. 在MATLAB 工作窗口输入程序

>>

[k,x,wuca,yx]=erfen (-2,-1,0.001)

运行后屏幕显示用二分法计算过程被列入表 2-3，其余结果为

k = 9，x=-1.7959，

wuca = 9.7656e-004，yx = 0.0037

图2-7

2.4 迭代法及其MATLAB 程序

确定根的近似位置以后，接下来的工作就是将根精确化，即按某种方法将初始近似值逐步精确化，直到满足所要求的精确度为止. 上节介绍的二分法是将根精确化的方法之一，但是它的收敛速度较慢，且不能求出偶重根.迭代法可以克服这种缺陷.迭代法是求解方程的根、线性和非线性方程组的解的基本而重要的方法.

2.4.2 迭代法的MATLAB 程序1

迭代法需要自行编制程序.下面提供的迭代法的MATLAB 程序1使用时只需输入迭代初始值0x 、迭代次数k 、迭代公式x k+1=?(x k )和一条命令，运行后就可以输出求迭代序列}{k x 、迭代k 次得到的迭代值k x 和相邻两次迭代的偏差piancha =|1--k k x x | （简称偏差）和偏差的相对误差xdpiancha =|1--k k x x |

k x 的值，并且具有警报功能（若迭代序列发散，

则提示用户“请重新输入新的迭代公式”；若迭代序列收敛，则屏幕会出现“祝贺您！此迭代序列收敛，且收敛速度较快”）.我们可以用这个程序来判断迭代序列的敛散性，也可以用于比较由一个方程得到的几个迭代公式的敛散性的优劣.

function [k,piancha,xdpiancha,xk]=diedai1(x0,k) % 输入的量--x0是初始值,k 是迭代次数 x(1)=x0; for i=1:k

x(i+1)=fun1(x(i));%程序中调用的fun1.m 为函数y=φ(x)

piancha= abs(x(i+1)-x(i)); xdpiancha= piancha/( abs(x(i+1))+eps); i=i+1;xk=x(i);[(i-1) piancha xdpiancha xk] end

if (piancha >1)&(xdpiancha>0.5)&(k>3)

disp('请用户注意：此迭代序列发散，请重新输入新的迭代公式') return ; end

if (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝贺您！此迭代序列收敛，且收敛速度较快') return ; end

p=[(i-1) piancha xdpiancha xk]';

例2.4.1 求方程102)(2-+=x x x f 的一个正根. 解 在MATLAB 工作窗口输入程序

>> [k,piancha,xdpiancha,xk]= diedai1(2,5)

运行后输出用迭代公式2/)10(2

1k k x x -=+的结果

[k,piancha,xdpiancha,xk]=

1.00000000000000 1.00000000000000 0.33333333333333 3.00000000000000

2.00000000000000 2.50000000000000 5.00000000000000 0.50000000000000

3.00000000000000

4.37500000000000 0.89743589743590 4.87500000000000

4.00000000000000 11.75781250000000 1.70828603859251 -6.88281250000000

5.00000000000000 11.80374145507813 0.63167031671297 -18.68655395507813 请用户注意：此迭代序列发散，请重新输入新的迭代公式

k=5，piancha = 11.80374145507813，xdpiancha = 0.63167031671297， xk = -18.68655395507813

由以上运行后输出的迭代序列与根*

x =2.316 624 790 355 40相差越来越大，即迭代序列}{k x 发散，此迭代法就失败.这时偏差piancha 逐渐增大且偏差的相对误差xdpiancha 的值大于0.5.

用迭代公式)2/(101+=+k k x x 运行后输出的结果 [k,piancha,xdpiancha,xk]=

1.00000000000000 0.50000000000000 0.20000000000000

2.50000000000000 2.00000000000000 0.27777777777778 0.12500000000000 2.22222222222222

3.00000000000000 0.14619883040936 0.06172839506173 2.36842105263158

4.00000000000000 0.07926442612555 0.03462603878116 2.28915662650602

5.00000000000000 0.04230404765128 0.01814486863115 2.33146067415730 k =5,piancha =0.04230404765128,xdpiancha = 0.01814486863115, xk = 2.33146067415730

可见，偏差piancha 和偏差的偏差的相对误差xdpiancha 的值逐渐变小，且第5次的迭代值x k =2.331 460 674 157 30与根*

x =2.316 624 790 355 40接近，则迭代序列}{k x 收敛，但收敛速度较慢，此迭代法较为成功.

用迭代公式)22/()102(2

1+-+-=+k k k k k x x x x x 运行后输出的结果

[k,piancha,xdpiancha,xk ]=

1.00000000000000 0.33333333333333 0.14285714285714

2.33333333333333

2.00000000000000 0.01666666666667 0.00719424460432 2.31666666666667

3.00000000000000 0.00004187604690 0.00001807631822 2.31662479061977

4.00000000000000 0.00000000026437 0.00000000011412 2.31662479035540

5.00000000000000 0 0 2.31662479035540 祝贺您！此迭代序列收敛，且收敛速度较快

k = 5; piancha = 0; xdpiancha = 0; y = 2.31662479035540.

可见，偏差piancha 和偏差的相对误差xdpiancha 的值越来越小，且第5次的迭代值x k =2.316 624 790 355 40与根*

x =2.316 624 790 355 40相差无几，则迭代序列}{k x 收敛，且收敛速度很快，此迭代法成功.

2.4.5 迭代法的MATLAB 程序2

function [k,piancha,xdpiancha,xk,yk]=diedai2(x0,tol,ddmax) x(1)=x0;

for i=1: ddmax

x(i+1)=fun(x(i));piancha=abs(x(i+1)-x(i)); xdpiancha=piancha/( abs(x(i+1))+eps);i=i+1;

xk=x(i);yk=fun(x(i)); [(i-1) piancha xdpiancha xk yk] if (piancha

if i>ddmax

disp('迭代次数超过给定的最大值ddmax') return ; end

P=[(i-1),piancha,xdpiancha,xk,yk]';

例2.4.5 求x 5-3x +1=0在0.3附近的根，精确到4位小数. 解 ⑴ 构造迭代公式)(k k x x ?=.

改写原方程为等价方程3/)1(5+=x x .这时3/)1()(5+=x x ?，初始值为x 0=0.5, 迭代公式

3/)1(5

1+=+k k x x ),2,1,0( =k .

⑵利用迭代法的MATLAB 程序2计算精确到4位小数的根的近似值 y 和迭代次数k .

在MATLAB 工作窗口输入程序

>> [k,piancha,xdpiancha,xk,yk]=diedai2(0.3,1e-4,100)

运行后输出的结果

[k,piancha,xdpiancha,xk,yk] =

0 0.03414 0.10218 0.30000 0.33414 1 0.03414 0.10218 0.33414 0.33472 2 0.00058 0.00173 0.33472 0.33473 3 0.00001 0.00004 0.33473 0.33473

k =3;piancha =1.206089525390697e-005; xdpiancha =3.603129477781680e-005;

xk =0.3347; yk =0.3347.

2.5 迭代过程的加速方法及其MATLAB 程序

对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢，从而使计算量变得很大，因此，迭代过程的加速是个重要的问题.

2.5.2 加权迭代法的MATLAB 程序

function [k,xk,yk]=jasudd (x0,tol,L,ddmax) x1(1)=x0;

for i=1: ddmax

x(i+1)=fun(x1(i));

x1(i+1)= x(i+1)+L*( x(i+1)- x1(i))/(1-L); piancha=abs(x1(i+1)-x1(i));

xdpiancha= piancha/( abs(x1(i+1))+eps); i=i+1;xk=x1(i);yk=fun(x1(i));[(i-1) xk yk] if (piancha

end end

if i>ddmax

disp('迭代次数超过给定的最大值ddmax') return ; end

P=[(i-1),xk,yk]';

例2.5.1 用（2.10）式求2e 0=--x x

在0.85附近的一个近似根，要求精度610-=ε. 解 在MATLAB 工作窗口输入程序

>> [k,xk,yk]=jasudd (0.85,1e-6,-0.855,100)

运行后输出结果

[k,xk,yk] =

1.00000000000000 0.85260370028041 0.85260703830561

2.00000000000000 0.85260549975491 0.85260550406236

3.00000000000000 0.85260550207699 0.85260550208255 k =3； xk =0.852606； yk = 0.852606.

2.5.4 艾特肯(Aitken)加速方法的MATLAB 程序

function [k,xk,yk,p]= Aitken (x0,tol, ddmax) x(1)=x0;

for i=1: ddmax

x1(i+1)=fun(x(i)); x2(i+1)=fun(x1(i+1));

x(i+1)=x2(i+1)-(x2(i+1)-x1(i+1))^2/(x2(i+1)-2*x1(i+1)+ x(i)); piancha=abs(x(i+1)-x(i));

xdpiancha= piancha/( abs(x(i+1))+eps); i=i+1; xk=x(i);yk=fun(x(i));

if (piancha

end

end

if i>ddmax

disp('迭代次数超过给定的最大值ddmax') return ; end

m=[0,1:i-1]; p=[m',x1',x2',x'];

例 2.5.3 用艾特肯加速方法求2e

0=--x x

在0.85附近的一个近似根，要求精度

610-=ε.

解 在MATLAB 工作窗口输入程序

>> [k,xk,yk,p]= Aitken(0.85,1e-6, 100)

运行后输出结果

k=3，xk=0.85260550201343，yk=0.85260550201373

p = 0 0 0 0.85000000000000 1.00000000000000 0.85482986389745 0.85071110652484 0.85260683568607 2.00000000000000 0.85260436491811 0.85260647150826 0.85260550201407 3.00000000000000 0.85260550201343 0.85260550201398 0.85260550201373

2.6 牛顿(Newton)切线法及其MATLAB 程序

2.6.2 牛顿切线法的收敛性及其MATLAB 程序

的供程序：

function [y,f]=newjushou(x)

f=fnq(x); fz=fnq(x)*ddfnq(x)/((dfnq(x))^2+eps); y=abs(fz); if (y<1)

disp('恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|d φ(x)/dx|和方程

f(x)=0的函数f(x)值f 如下')

else

disp('请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程

f(x)=0的函数f(x)值')

end

P=[y,f]';

例 2.6.2 用牛顿切线法的局部收敛性判别方程 e 4sin =x x 的近似根时，由下列初始值0x 产生的迭代序列是否收敛？

⑴-1;0=x ⑵0;0=x ⑶1;0=x ⑷2;0=x ⑸ 5.50=x ；⑹80=x . 解 在MATLAB 工作窗口输入程序

>> [y,f]=newjushou(-1)

运行后输出结果

请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值

y =

139.5644 f =

4.3096

说明：实际上，这个初始值所产生的迭代序列{}k x 发散,见例 2.6.6. （2）输入程序

>> [y,f]=newjushou(0)

运行后输出结果

请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值

y =

8.0000 f = 4

说明：实际上，这个初始值所产生的迭代序列{}k x 收敛到方程的根2.925 32，而不是收敛到离它最近的根1.400 81,即不是局部收敛的，见例 2.6.6.

（3）输入程序

>> [y,f]=newjushou(1)

运行后输出结果

恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值f 如下

y =

0.3566

f =

1.7126

说明：实际上，这个初始值所产生的迭代序列{}k x 收敛到离它最近的根1.400 81,即{}k x 是局部收敛的，见例 2.6.6. （4）输入程序：

>> [y,f]=newjushou(2)

运行后输出结果

请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值

y =

1.2593 f =

-2.7188

说明：虽然y =1.259 31>，不满足牛顿切线法局部收敛的条件1

)(*

<'x ? .但是，实

际上，这个初始值所产生的迭代序列{}k x 收敛到离它最近的根1.400 81,即{}k x 是局部收敛的.这说明1)(0<'x ?是充分条件，非必要条件（见例 2.6.6）. （5）输入程序

>> [y,f]=newjushou(5.5)

运行后输出结果

请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值

y =

1.0447e+005 f =

176.6400

说明：实际上，这个初始值所产生的迭代序列{}k x 收敛到235.619,但不是此方程的根，

比较初始值5.5与迭代值235.619，可见{}k x 不是局部收敛的，见例 2.6.6. （6）输入程序

>> [y,f]=newjushou(8)

运行后输出结果

恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值f 如下

y =

0.4038 f =

-2.9452e+003

说明：实际上，这个初始值所产生的迭代序列{}k x 收敛到方程的根6.290 60，而不是收敛到离它最近的根9.424 46,即不是局部收敛的，见例 2.6.6.

2.6.3 牛顿切线法的MATLAB 程序

牛顿切线法求方程0)(=x f 的近似根x （精度为tol ）需要自行编制程序.

function [k,xk,yk,piancha,xdpiancha]=newtonqx(x0,tol,ftol,gxmax)

x(1)=x0;

for i=1: gxmax

x(i+1)=x(i)-fnq(x(i))/(dfnq(x(i))+eps); piancha=abs(x(i+1)-x(i));

xdpiancha= piancha/( abs(x(i+1))+eps); i=i+1;

xk=x(i);yk=fnq(x(i)); [(i-1) xk yk piancha xdpiancha] if (abs(yk)

if i>gxmax

disp('请注意：迭代次数超过给定的最大值gxmax 。') return ; end

[(i-1),xk,yk,piancha,xdpiancha]';

例 2.6.3 用牛顿切线法求方程01322

3=+-x x 在9.04.00和-=x 的近似根，要求精

度3

10-=ε.然后判断此方法分别在方程的根15.0*

和-=x 处的收敛速度.

解 （1）先求方程的近似根. 在MATLAB 工作窗口输入程序

>> [k,xk,yk,piancha,xdpiancha]=newtonqx(-0.4,0.001, 0.001,100)

运行后输出初始值4.00-=x 的迭代结果如表 2-10所示.

k k =1.712e-007和相对误差xdpiancha =3.423e-007.

如果将步骤④命令中的-0.4改作0.9，则运行后输出初始值4.00-=x 的迭代结果为 k =7；piancha =7.290e-004；xdpiancha =7.295e-004；xk = 0.999；yk =1.590e-006.

（2）再讨论收敛速度

比较初始值分别是9.04.00和-=x 的两组结果，在4.00-=x 处迭代3次就得到单根

5.0*-=x 的近似值x k =-0.500；而在9.00=x 处迭代7次才得到二重根1*=x 的近似值

x k =0.999.可见，牛顿切线迭代法在单根处的迭代速度比二重根处的迭代速度快很多.这正如前面讨论的结果一样，即若*

x 是0)(=x f 的单根，则牛顿切线法是平方收敛的；若*

x 是

0)(=x f 的二重根，

则牛顿切线法是线性收敛的；我们也可以用阶的定义判断它们的收敛速度.留给读者证明.

2.6.4 求n c 的方法及其MATLAB 程序

用c 的n 次根n c 的迭代公式求n c 的近似值x （精度为tol ）需要自行编制程序.

根据求c 的n 次根n c 的迭代公式(2.25)和已知条件，现提供名为kai2fang .m 的M 文件： function [k,xk,yk,piancha,xdpiancha,P]=kainfang(x0,c,n,tol, gxmax)

x(1)=x0; for i=1: gxmax

u(i)= (x(i)^n-c)/(n*x(i)^(n-1)); x(i+1)= x(i)-u(i); piancha=abs(x(i+1)-x(i)); xdpiancha=piancha/( abs(x(i+1))+eps); i=i+1; xk=x(i);yk=fnq(x(i));

[(i-1),xk,yk,piancha,xdpiancha]

if (piancha

end end

if i>gxmax

disp('请注意：迭代次数超过给定的最大值gxmax.')

[(i-1),xk,yk,piancha,xdpiancha] return ; end

P=[(i-1),xk,yk,piancha,xdpiancha]';

例2.6.5 求113，要求精度为6

10-.

解 本题介绍四种解法.

方法1 用求c 的n 次根n c （当n 是偶数时，0>c ）的MATLAB 程序计算.

取初始值100=x ，根指数2=n ，被开方数113=c ，近似根的精度tol =5

10-，迭代

的最大次数gxmax =100.在工作区间输入程序

>> [k,xk,yk,piancha,xdpiancha]=kainfang(10,113,2,1e-5,100)

运行后输出结果

k=4,piancha=1.610800381968147e-011, xdpiancha=1.515313534117706e-012

xk =10.63014581273465,yk =1.710910332060509e+009 可见,≈11310.630 15，满足精度5

10-.

方法2 用牛顿迭代公式（2.12）计算.

设113=x ，则01132

=-x ，记113)(2-=x x f ，x x f 2)('=.由牛顿迭代公式得，

k k k k x x x x 21132

1--

=+，即）（k

k k x x x 113

211+=+)321,0( ，，，其中=k 取初始值100=x ，计算结果列入表 2-12.

因为，迭代次数k =4时，偏差34x x -810-<，满足精度5

10-，所以，

≈11310.630 15.

方法3 用牛顿切线法的MATLAB 主程序计

算.

分别建立名为fnq.m 和dfnq.m 的M 文件

function y=fnq(x)

y=x^2-113;

function y=dfnq(x)

y=2*x;

在MATLAB 工作窗口输入程序

>> [k,xk,yk,piancha,xdpiancha]=newtonqx(10, 1e-5, 1e-5,100)

运行后，将输出的结果列入下表 2-13. 迭代k =4次，得到精度为5

10-的结果≈113

10.630 15.

方法4 在MATLAB 工作空间输入程序

>> 113^(1/2)

运行后输出

ans = 10.63014581273465 经过四舍五入后，得到精度为5

10-的结果≈11310.630 15.

2.6.6 牛顿切线法的加速及其两种MATLAB 程序

当*x 是方程0)(=x f 的)2(≥m m 重根时，由牛顿切线公式（2.12）产生的迭代序列{}k x 收敛到*

x 的速度较慢，是线性收敛的.我们要提高收敛速度，可以有选择地采用已

知根的重数和未知根的重数两种求重根的修正牛顿切线公式.

function [k,piancha,xdpiancha,xk,yk]=newtonxz(m,x0,tol,ftol,gxmax) x(1)=x0;

for i=1: gxmax

x(i+1)=x(i)-m*fnq(x(i))/(dfnq(x(i))+eps); piancha=abs(x(i+1)-x(i));

xdpiancha=piancha/( abs(x(i+1))+eps); i=i+1; xk=x(i);yk=fnq(x(i));

if ((piancha

k=i-1; xk=x(i);[(i-1) piancha xdpiancha xk yk] return ; end end

if i>gxmax

disp('请注意：迭代次数超过给定的最大值gxmax.') return ; end

P=[(i-1),piancha,xdpiancha,xk,yk]';

例 2.6.7 判断0*=x 是方程2)(=x f e 02692

3=---x x x 的几重根？在区间

]1,0[上，分别用牛顿切线法和求重根的修正牛顿切线公式求此根的近似值k x ，使精确到

410-=ε.

解 经判断知，0*=x 是方程0)(=x f 的三重根.在MATLAB 工作窗口输入程序

>> [k,piancha,xdpiancha,xk,yk]= newtonqx (0.5,1e-4, 1e-4,100)

运行后整理结果列入表 2-20.

迭代次数k =20，精确到4

10-=ε的根的近似值是x k =0.000 2，其函数值是y k =5.755 8e-011，是一阶（线性）收敛速度.

根据求重根的修正牛顿切线公式，在MATLAB 工作窗口输入程序

>> [k,piancha,xdpiancha,xk,yk]= newtonxz(3,0.5,1e-4, 1e-4,100)

运行后整理结果得表 2-21.

迭代次数k =7，精确到4

10-=ε的根的近似值是x k =1.121 1e-006，piancha = 9.975

3e-005，xdpiancha =88.974 8，其函数值是y k = -8.881 8e-016，是二阶收敛（平方收敛）速度.可见，求重根的修正牛顿切线公式比牛顿切线法收敛速度快得多.

function [k,piancha,xdpiancha,xk,yk]=newtonxz1(x0,tol,ftol,gxmax) x(1)=x0;

for i=1: gxmax

u(i)=fnq(x(i))/dfnq(x(i));

du(i)=1-fnq(x(i))*ddfnq(x(i))/((dfnq(x(i)))^2+eps); x(i+1)=x(i)-u(i)/du(i); piancha=abs(x(i+1)-x(i));

xdpiancha=piancha/( abs(x(i+1))+eps); i=i+1; xk=x(i);yk=fnq(x(i)); if ((piancha

k=i-1; xk=x(i); [(i-1) piancha xdpiancha xk yk] return ; end end

if i>gxmax

disp('请注意：迭代次数超过给定的最大值gxmax.') return ; end

P=[(i-1),piancha,xdpiancha,xk,yk]';

例2.6.8 用未知重数的求重根的修正牛顿切线法，求方程 在区间]1,0[上的根，使精确到4

10-=ε，并且与例2.6.7的结果比较. 解 在MATLAB 工作窗口输入程序

>> [k,piancha,xdpiancha,xk,yk]=newtonxz1(0.5,1e-4, 1e-4,100) 运行后整理结果得表 2-22.

即迭代k =4次，piancha =2.461 55e-005， xdpiancha =0.945 62， x k =-2.603 09e-005， y k =-1.325 61e-013.这些结果与例2.6.7的结果比较见下表 2-23.

可见，未知重数的修正牛顿切线法的迭代次数和偏差最小；已知重数的修正牛顿切线法的迭代次数比牛顿切线法少13次，比未知重数的修正牛顿切线法多3次，但是根的近似值最接近根0；用牛顿切线法求得的近似根是三者中与根0距离最远的.

2.7 割线法及其MATLAB 程序

割线法又称弦截法、弦位法、弦割法和弦法.

牛顿切线法的突出优点是收敛速度快，但它也有明显的缺点：公式中含有导数，当)(x f 较复杂时，使用不方便.为了避免切线法计算导数的麻烦，我们现在设法利用一些函数值)(k x f ，)(1-k x f ，…来回避导数值)('k x f 的计算.下面具体介绍割线法.

2.7.2 割线法的MATLAB 程序

用割线法求方程0)(=x f 的近似根x （精度为tol ）需要自行编制程序.

使用初始值0和1，根据式（2.28），现提供名为gexian.m 的M 文件：

function [k,piancha,xdpiancha,xk,yk]=gexian (x01,x02,tol,ftol,gxmax) x(1)=x01;x(2)=x02; for i=2: gxmax

u(i)= fnq(x(i))*(x(i)-x(i-1)); v(i)= fnq(x(i))-fnq(x(i-1)); x(i+1)=x(i)- u(i)/( v(i)); piancha=abs(x(i+1)-x(i));

xdpiancha= piancha/( abs(x(i+1))+eps); i=i+1; xk= x(i);

yk=fnq(x(i)); [(i-2) piancha xdpiancha xk yk] if (abs(yk)

if i>gxmax

disp('请注意：迭代次数超过给定的最大值gxmax.') return ; end

P=[(i-2),piancha,xdpiancha,xk,yk]';

例2.7.1 用割线法求方程=)(x f e

0322=-x x

的一个实根的近似值k x ，使精确到

410-=ε.

解 先用作图法确定初始值x 01和x 02，然后在MATLAB 工作窗口输入程序

>> [k,piancha,xdpiancha,xk,yk]= gexian (-0.4,-0.3,1e-4, 1e-4,100)

运行后输出结果经整理得表 2-24.由此可见，迭代k =3次，得到精度为3

10-的根的近似值x k =-0.390 6，其函数值为y k =3.860 7e-008，x k 的相邻偏差为piancha =3.327 1e-005，其相误差为xdpiancha =8.516 8e-005.

2.8 抛物线法及其MATLAB 程序

2.8.2 抛物线法的MATLAB 程序

用抛物线法求方程0)(=x f 的近似根

x （精度为tol ）需要自行编制程序.

使用初始值0，1和2，根据式按公式（2.35），（2.36）和（2.39），现提供名为paowu.m 的M 文件

function [k,piancha,xdpiancha,xk,yk]=paowu(x1,x2, x3,tol,ftol,gxmax)

X=[x1, x2, x3]; for i=1: gxmax

h0= X(1)- X (3); h1= X (2)- X (3); u0= fnq (X (1))- fnq (X (3) ); u1= fnq (X (2))- fnq (X (3)); c= fnq (X (3));

fenmu= h0^2* h1- h1^2* h0; afenzi= u0* h1- u1* h0; bfenzi= u1*h0^2-u0*h1^2;

a= afenzi/ fenmu; b= bfenzi/ fenmu; deta=b^2-4*a*c;

cha=[abs(X(3)-X(1)),abs(X(3)-X(2)),abs(X(2)- X(1))];

piancha =min(cha); xdpiancha= piancha/( abs (X(3))+eps); if deta<0

disp('提示：根的判别式值为负数.'),detakf=sqrt(deta); elseif deta<=0

disp('提示：根的判别式值为零.'),detakf=0; else

disp('提示：根的判别式值为正数.'),detakf=sqrt(deta); end if b<0

detakf=- detakf; end

z=-2*c/(b+ detakf);

xk= X(3)+ z;q1= xk - X (1); q2= xk - X (2); q3= xk - X (3); if abs(q2)< abs(q1)

X1=[X(2), X(1), X(3)]; X= X1;Y= fnq(X); end

if abs(q3)< abs(q1)

X2=[X(1), X(3), X(2)]; X= X2;Y= fnq(X);

end

X(3)= xk; Y(3)=fnq(X(3));

yk= Y(3); [i piancha xdpiancha xk yk]

if (abs(yk)

if i>gxmax

disp('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值x0') return ; end

P=[i,piancha,xdpiancha,xk,yk]';

例2.8.1 用抛物线法求方程=)(x f e

0322=-x x

的一个实根的近似值k x ，使精确到

410-=ε.

解 先用作图法确定初始值x 01，x 02和x 03，然后在MATLAB 工作窗口输入程序

>> [k,piancha,xdpiancha,xk,yk]= paowu (-0.4,-0.3, -0.5,1e-4, 1e-4,100)

运行后输出结果

提示：根的判别式值为正数. ans =

Columns 1 through 4

1.00000000000000 0.10000000000000 0.20000000000000 -0.39066350562239 Column 5

-0.00005581900299

提示：根的判别式值为正数. ans =

Columns 1 through 4

2.00000000000000 0.00933649437761 0.02389906977038 -0.39064638365394 Column 5

-0.00000000923532

提示：根的判别式值为正数. ans =

Columns 1 through 4

3.00000000000000 0.00001712196845 0.00004382983990 -0.39064638082059 Column 5

0.00000000000000 k = 3 piancha =

1.712196845282676e-005 xdpiancha =

4.382983989938760e-005 xk =

-0.39064638082059 yk =

2.220446049250313e-016

即，迭代k =3次，得到精度为4

10-的根的近似值x k =-0.390 6，其函数值为y k =2.220 4e-016，x k 的相邻最小偏差为piancha =1.712e-005，其相对误差为xdpiancha =4.383 0e-005.

将抛物线法与割线法的迭代结果进行比较，见表 2-25.显然，抛物线法比割线法的迭代速度快.

2.9 求解非线性方程组的牛顿法及其MATLAB 程序

求解非线性方程的牛顿法可以推广到二维非线性方程组的情形，也可以推广到更高维的情形.

2.9.1 求解二元非线性方程组的牛顿法及其MATLAB 程序

用牛顿法求二元非线性方程组（2.40）的近似根k k y x ,（精度为tol ）需要自行编制程

使用两个初始值00，根据公式（2.45），（2.46），（2.47）和（2.49），现提供名为newtonzu2.m 的M 文件

function [ci,D,danfan,xddf,hanfan,Xk,Yk]=newtonzu2(X,tol,ftol,ngmax) Y=Z(X);

for i=1:ngmax

dY=JZ(X);D=det(dY);

A1=[Y(1),dY(1,2); Y(2),dY(2,2)]; A=det(A1);

B1=[dY(1,1), Y(1); dY(2,1), Y(2)];B=det(B1);Xk=X-[A,B]/D; hanfan=norm(Y);danfan=norm(Xk-X); xddf=danfan/(norm(Xk)+eps); X=Xk;Y=Z(X);ci=i; if D~=0

ci=i; Xk=X-[A,B]/D;

Yk=Y; [ci,D,danfan, xddf, hanfan , X, Y]; else

disp('请注意！迭代方程组（2.44）的系数行列式的值等于零.') end

if (hanfan

[ci,D,danfan, xddf, hanfan , X, Y]; return ; end

if i>gxmax

disp('请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值x0.') return ; end end

例2.9.1 用迭代公式（2.46）求解方程组 ???

??=-=+.

2,162222y x y x

解 这里介绍两种求解已知方程组的方法.

方法1 用解二元非线性方程组的牛顿法的MATLAB 主程序求解.

（1）设16)(2

21-+=y x x f ，

=)(2x f 222--y x ，

则

x x

y x f 2)

,(1=??，

y y

y x f 2)

,(1=??，

x x y x f 2)

,(2=??，

y y

y x f 2)

,(2-=??.

（2）作函数)(

1x f 和)(2x f 的图形.输入程

序

>>syms x y

F1=x^2-y^2-2;F2=x^2+y^2-16; ezplot(F1,[-4.5,4.5]), hold on ezplot(F2,[-4.5,4.5]) ,hold off

运行后屏幕显示图2-29.取两个初始值,20=x 20=y .

（3）建立并保存名为Z.m 的M 文件

function F=Z(X)

x=X(1);y=X(2); F=zeros(1,2);

F(1)= x^2+y^2-16; F(2)= x^2-y^2-2; （4）建立并保存名为JZ.m 的M 文件

function dF=JZ(X)

x=X(1);y=X(2); dF=zeros(2,2); dF(1,1)=2*x; dF(1,2)=2*y; dF(2,1)=2*x; dF(2,2)=-2*y;

（5）保存名为newtonzu2.m 的M 文件； 图2-29 （6）在MATLAB 工作窗口输入程序:

>> X=[2,2]; [k,D,danfan, xddf, hanfan , Xk, Yk]= newtonzu2 (X,1e-4,1e-4,100) （7）运行后输出结果

k = 5 D =

-63.49803146638493 danfan =

3.932201289951168e-011 xddf =

9.830503224877920e-012 hanfan =

3.336593498083849e-010 Xk =

2.99999999996068 2.64575131106449 Yk =

1.0e-015 *

0 -0.88817841970013 方法2 用解矩阵方程的方法求解. 编写如下程序

>> x0=2;y0=2; gxmax =10;tol=1e-6;x(1)=x0;y(1)=y0; i=1;u=[1 1];k(1)=1; while (norm(u)>tol*norm([x(i),y(i)]’)) A=2*[x(i),y(i);x(i),-y(i)];

b=[16-x(i)^2-y(i)^2,2-x(i)^2+y(i)^2]’;

u=A\b;x(i+1)=x(i)+u(1); y(i+1)=y(i)+u(2); i=i+1;k(i)=i; if (i> gxmax)

error('请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值'); end end

[k’,x’,y’] （2）运行后输出结果,取初始值=x T

)2,2(，迭代次数n =6，精度ε=10-6

，