本章主要介绍方程根的有关概念

1、第二章 非线性方程（组）的数值解法本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（Müller）法和迭代法的加速等）及其MATLAB程序，求解非线性方程组的方法及其MATLAB程序.2.1 方程（组）的根及其MATLAB命令求解方程（组）的solve命令求方程f(x)=q(x)的根可以用MATLAB命令:>> x=solve('方程f(x)=q(x)',待求符号变量x)求方程组fi(x1,xn)=qi(x1,xn) (i=1,2,n)的根可以用MATLAB命令

2、:>>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)求解多项式方程（组）的roots命令如果为多项式，则可分别用如下命令求方程的根，或求导数（见表 2-1）.表 2-1 求解多项式方程（组）的roots命令命 令功 能xk =roots(fa)输入多项式的系数fa（按降幂排列），运行后输出xk为的全部根.dfa=polyder(fa)输入多项式的系数fa（按降幂排列），运行后输出dfa为多项式的导数的系数.

3、dfx=poly2sym(dfa)输入多项式的导数的系数dfa（按降幂排列），运行后输出dfx为多项式的导数. 求解方程（组）的fsolve命令如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots命令.如果非线性方程（组）是含有超越函数，则无法使用roots命令，需要调用MATLAB系统中提供的另一个程序fsolve来求解.当然，程序fsolve也可以用于多项式方程（组），但是它的计算量明显比roots命令的大. fsolve命令使用最小二乘法（least squares method）解非线性方程（组）的数值解，其中X和F（X）可以是向量或矩阵.此

4、种方法需要尝试着输入解X的初始值（向量或矩阵）X0，即使程序中的迭代序列收敛，也不一定收敛到的根（见例2.1.8）.fsolve的调用格式： X=fsolve(F,X0)输入函数的M文件名和解X的初始值（向量或矩阵）X0，尝试着解方程（组），运行后输出解的估计值（向量或矩阵）X.要了解更多的调用格式和功能请输入:help fsolve,查看说明.2.2 搜索根的方法及其MATLAB程序求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根，需要搜索根所在的区间或确定根的初始近似值（简称初始值）.搜索根的近似位置的常用方法有三种:作图法、逐步搜索法和二分法等，使用这些

5、方法的前提是高等数学中的零点定理. 作图法及其MATLAB程序作函数的图形的方法很多，如用计算机软件的图形功能画图，或用高等数学中应用导数作图，或用初等数学的函数叠加法作图等.下面介绍两种作图程序.作函数在区间 a,b 的图形的MATLAB程序一x=a:h:b; % h是步长y=f(x); plot(x,y)grid, gtext('y=f(x)')说明： 此程序在MATLAB的工作区输入，运行后即可出现函数的图形.此图形与轴交点的横坐标即为所要求的根的近似值. 区间a,b 的两个端点的距离 b-a 和步长h的绝对值越小，图形越精确.作函数在区间 a,b上的图形的MATLAB程

6、序二将化为，其中是两个相等的简单函数x=a:h:b; y1=h(x); y2=g(x);plot(x, y1, x, y2)grid,gtext(' y1=h(x),y2=g(x)')说明：此程序在MATLAB的工作区输入，运行后即可出现函数的图形.两图形交点的横坐标即为所要求的根的近似值.下面举例说明如何用计算机软件MATLAB的图形功能作图. 逐步搜索法及其MATLAB程序逐步搜索法也称试算法.它是求方程根的近似值位置的一种常用的方法. 逐步搜索法依赖于寻找连续函数满足与异号的区间.一旦找到区间，无论区间多大，通过某种方法总会找到一个根.MATLAB的库函数中没有逐步搜索法

7、的程序，现提供根据逐步搜索法的计算步骤和它的收敛判定准则编写其主程序，命名为zhubuss.m.逐步搜索法的MATLAB主程序 输入区间端点a和b的值，步长和精度tol，运行后输出迭代次数 k=(b-a)/h+1,方程 根的近似值r.function k,r=zhubuss(a,b,h,tol)% 输入的量- 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)

8、;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)<tol)&( xielv<=0) m=m+1;r(m)=X(k); endend例 用逐步搜索法的MATLAB程序分别求方程和在区间上的根的近似值，要求精度是0.000 1.解 在MATLAB工作窗口输入如下程序>> k,r=zhubuss(-2,2,0.001,0.

9、0001)运行后输出的结果k =4001r = -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) 代替，可得到方程在区间上的根的近似值如下r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.9230 1.3310 1.5780 1.7650 1.9200如果读

10、者分别将方程的结果与例2.2.3比较，方程的结果与例比较，将会发现逐步搜索法的MATLAB程序的优点.如果精度要求比较高，用这种逐步搜索法是不合算的.2.3 二分法及其MATLAB程序 二分法二分法也称逐次分半法.它的基本思想是：先确定方程含根的区间（a,b），再把区间逐次二等分.我们可以根据式（2.3b）、区间a,b和误差，编写二分法求方程根的迭代次数的MATLAB命令.已知闭区间a,b和误差.用二分法求方程误差不大于的根的迭代次数的MATLAB命令k=-1+ceil(log(b-a)- log(abtol)/ log(2) % ceil是向+方向取整，abtol是误差. 二分法的MATLA

11、B程序二分法需自行编制程序，现提供用二分法求方程f(x)=0的根的近似值的步骤和式（2.3a）编写一个名为erfen.m的二分法的MATLAB主程序如下.二分法的MATLAB主程序 求解方程在开区间（a,b）内的一个根的前提条件是在闭区间a,b上连续，且.输入的量：a和b是闭区间a,b的左、右端点，abtol是预先给定的精度.运行后输出的量：k是使用二分法的次数.x是方程 在（a,b）内的实根x*的近似值，其精度是abtol.wuca=|bk-ak|/2是使用k次二分法所得到的小区间的长度的一半，即实根x*的近似值x的绝对精度限,满足wuca abtol.yx=f(xk) ，即方程f(x)=0

12、在实根x*的近似值x处的函数值.function k,x,wuca,yx=erfen(a,b,abtol)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.'), returnendmax1=-1+ceil(log(b-a)- log(abtol)/ log(2); % ceil是向 方向取整for k=1: max1+1a;ya=fun(a); b;yb=fun(b); x=(a+b)/2; yx=fun(x);

13、 wuca=abs(b-a)/2; k=k-1;k,a,b,x,wuca,ya,yb,yxif yx=0a=x; b=x;elseif yb*yx>0b=x;yb=yx;elsea=x; ya=yx;endif b-a< abtol , return, endendk=max1; x; wuca; yx=fun(x);例2.3.3 确定方程x3-x+4=0的实根的分布情况，并用二分法求在开区间 (-2,-1)内的实根的近似值，要求精度为0.001.图2-7解 （1）先用两种方法确定方程x3-x+4=0 的实根的分布情况。方法1 作图法.在MATLAB工作窗口输入如下程序>&g

14、t;x=-4:0.1:4; y=x.3-x +4; plot(x,y) grid,gtext('y=x3-x+4') 画出函数f(x)=x3-x+4的图像，如图2-7.从图像可以看出，此曲线有两个驻点都在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)满足，所以此方程在(-2,-1

15、)内有一个实根.（2） 用二分法的主程序计算.在MATLAB工作窗口输入程序>> k,x,wuca,yx=erfen (-2,-1,0.001)运行后屏幕显示用二分法计算过程被列入表 2-3，其余结果为k = 9，x=-1.7959， wuca = 9.7656e-004，yx = 0.00372.4 迭代法及其MATLAB程序确定根的近似位置以后，接下来的工作就是将根精确化，即按某种方法将初始近似值逐步精确化，直到满足所要求的精确度为止. 上节介绍的二分法是将根精确化的方法之一，但是它的迭代法的MATLAB程序1迭代法需要自行编制程序.下面提供的迭代法的MATLAB程序1使用时只

16、需输入迭代初始值、迭代次数k、迭代公式xk+1=(xk)和一条命令，运行后就可以输出求迭代序列、迭代k次得到的迭代值和相邻两次迭代的偏差piancha =| （简称偏差）和偏差的相对误差xdpiancha=|的值，并且具有警报功能（若迭代序列发散，则提示用户“请重新输入新的迭代公式”；若迭代序列收敛，则屏幕会出现“祝贺您！此迭代序列收敛，且收敛速度较快”）.我们可以用这个程序来判断迭代序列的敛散性，也可以用于比较由一个方程得到的几个迭代公式的敛散性的优劣.迭代法的MATLAB程序1 输入的量：初始值、迭代次数k和迭代公式；运行后输出的量：迭代序列、迭代k次得到的迭代值、相邻两次迭代的偏差pia

17、ncha =|和它的偏差的相对误差xdpiancha=的值.根据迭代公式（2.4）和已知条件，现提供名为diedai1.m的M文件如下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 xkendif

18、(piancha >1)&(xdpiancha>0.5)&(k>3) disp('请用户注意：此迭代序列发散，请重新输入新的迭代公式') return; end if (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝贺您！此迭代序列收敛，且收敛速度较快') return; endp=(i-1) piancha xdpiancha xk' 求方程的一个正根.解 在MATLAB工作窗口输入程序>> k,piancha,x

19、dpiancha,xk= diedai1(2,5)运行后输出用迭代公式的结果k,piancha,xdpiancha,xk=1.00000000000000 1.00000000000000 0.33333333333333 3.000000000000002.00000000000000 2.50000000000000 5.00000000000000 0.500000000000003.00000000000000 4.37500000000000 0.89743589743590 4.875000000000004.00000000000000 11.75781250000000 1.7

20、0828603859251 -6.88281250000000请用户注意：此迭代序列发散，请重新输入新的迭代公式k=5，piancha = 11.80374145507813，xk = -18.68655395507813由以上运行后输出的迭代序列与根=2.316 624 790 355 40相差越来越大，即迭代序列发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5.用迭代公式运行后输出的结果k,piancha,xdpiancha,xk= 1.00000000000000 0.50000000000000 0.20000000000000 2.

21、50000000000000k =5,piancha =0.04230404765128,可见，偏差piancha和偏差的偏差的相对误差xdpiancha的值逐渐变小，且第5次的迭代值xk =2.331 460 674 157 30与根=2.316 624 790 355 40接近，则迭代序列收敛，但收敛速度较慢，此迭代法较为成功.用迭代公式运行后输出的结果k,piancha,xdpiancha,xk= 祝贺您！此迭代序列收敛，且收敛速度较快.可见，偏差piancha和偏差的相对误差xdpiancha的值越来越小，且第5次的迭代值xk=2.316 624 790 355 40与根=2.316

22、624 790 355 40相差无几，则迭代序列收敛，且收敛速度很快，此迭代法成功.迭代法的MATLAB程序2迭代法的MATLAB程序2 输入的量：初初始值、精度tol和和迭代公式；运行后输出的量：迭代次数k，满足精度tol的根的近似根的迭代序列，相邻两次迭代的偏差pianch=|和它的偏差的相对误差xdpiancha =的值，yk=(xk).根据迭代公式（2.4）和已知条件，现提供名为diedai2.m的M文件如下function k,piancha,xdpiancha,xk,yk=diedai2(x0,tol,ddmax)x(1)=x0;for i=1: ddmax x(i+1)=fun(

23、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<tol)|(xdpiancha< tol) k=i-1; xk=x(i); return;endend if i>ddmax disp('迭代次数超过给定的最大值ddmax') return;endP=(i-1),piancha,xdpiancha,xk,yk' 求x5-3x+1=0在

24、0.3附近的根，精确到4位小数.解 构造迭代公式.改写原方程为等价方程.这时，初始值为x0=0.5, 迭代公式   .利用迭代法的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.334141 0.03414 0.10218 0.33414 0.334722 0.00058 0.00173 0

25、.33472 0.334733 0.00001 0.00004 0.33473 0.33473k =3;piancha =1.206089525390697e-005; xk =0.3347; yk =0.3347.2.5 迭代过程的加速方法及其MATLAB程序 对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢，从而使计算量变得很大，因此，迭代过程的加速是个重要的问题.MATLAB程序加权迭代法的MATLAB主程序已知初始值、精度tol和迭代公式xk+1=(xk)，求满足精度tol的根的近似根和它的函数值yk及迭代次数k.根据加权迭代的公式（2.10）和

26、已知条件，现提供名为jasudd.m的M文件如下：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<tol)|(xdpiancha< tol) k=i-1;

27、xk=x1(i); return;endendif i>ddmax disp('迭代次数超过给定的最大值ddmax') return;endP=(i-1),xk,yk' 用（2.10）式求2e在0.85附近的一个近似根，要求精度.解 在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.85

28、260550406236 3.00000000000000 0.85260550207699 0.85260550208255k =3； xk =0.852606； yk = 0.852606.MATLAB程序艾特肯加速方法的MATLAB程序已知初始值、精度tol和迭代公式xk+1=(xk)，求满足精度tol的根的近似根和它的函数值yk及迭代次数k，p=k',x1',x2',x'.根据艾特肯加速方法的公式（2.11）和已知条件，现提供名为Aitken.m的M文件如下：function k,xk,yk,p= Aitken (x0,tol, ddmax)x(1)=x

29、0;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<tol)|(xdpiancha<tol) k=i-1; xk=x1(i); return;endendif i>ddmax disp('迭

30、代次数超过给定的最大值ddmax') return;endm=0,1:i-1; p=m',x1',x2',x' 用艾特肯加速方法求2e在0.85附近的一个近似根，要求精度.解 在MATLAB工作窗口输入程序>> k,xk,yk,p= Aitken(0.85,1e-6, 100)运行后输出结果k=3，xk=0.85260550201343，yk=0.85260550201373p = 0 0 0 0.85000000000000 1.00000000000000 0.85482986389745 0.85071110652484 0.85260

31、683568607 2.00000000000000 0.85260436491811 0.85260647150826 0.85260550201407 3.00000000000000 0.85260550201343 0.85260550201398 0.852605502013732.6 牛顿(Newton)切线法及其MATLAB程序牛顿切线法的收敛性及其MATLAB程序判别牛顿切线法的局部收敛性的MATLAB主程序根据牛顿切线法局部收敛的条件 输入的量：x是方程根或满足的初始值，函数及其一阶导数和二阶导数；输出的量：或的值.如果或时，运行后输出信息'恭喜您！此迭代序列收敛,(

32、x)的导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值f如下'；否则，运行后输出信息'请注意观察下面显示的 (x)的导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值'.根据牛顿切线法局部收敛的条件和方程的函数及其一、二阶导数，现提供程序：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的函数

33、f(x)值f如下')else disp('请注意观察下面显示的(x)的导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值')end P=y,f' 用牛顿切线法的局部收敛性判别方程 e的近似根时，由下列初始值产生的迭代序列是否收敛？ ；.解 在MATLAB工作窗口输入程序>> y,f=newjushou(-1)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值y = 139.5644f = 4.3096说明：实际上，这个初始值所产生的迭代序列发散,见（2）输入程序 >

34、;> y,f=newjushou(0)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值y = 8.0000f = 4说明：实际上，这个初始值所产生的迭代序列收敛到方程的根2.925 32，而不是收敛到离它最近的根1.400 81,即不是局部收敛的，见（3）输入程序>> y,f=newjushou(1)运行后输出结果恭喜您！此迭代序列收敛,(x)导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值f如下y = 0.3566f = 1.7126说明：实际上，这个初始值所产生的迭代序列收敛到离它最近的

35、根1.400 81,即是局部收敛的，见（4）输入程序：>> y,f=newjushou(2)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值y = 1.2593f = -2.7188说明：虽然y =1.259 3，不满足牛顿切线法局部收敛的条件 .但是，实际上，这个初始值所产生的迭代序列收敛到离它最近的根1.400 81,即是局部收敛的.这说明是充分条件，非必要条件（见 （5）输入程序 >> y,f=newjushou(5.5)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d(x)/dx|和

36、方程f(x)=0的函数f(x)值y = 1.0447e+005f = 176.6400说明：实际上，这个初始值所产生的迭代序列收敛到235.619,但不是此方程的根，比较初始值5.5与迭代值235.619，可见不是局部收敛的，见（6）输入程序 >> y,f=newjushou(8)运行后输出结果恭喜您！此迭代序列收敛,(x)导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数f(x)值f如下y = 0.4038f = -2.9452e+003说明：实际上，这个初始值所产生的迭代序列收敛到方程的根6.290 60，而不是收敛到离它最近的根9.424 46,即不是局部收敛的，见

37、牛顿切线法的MATLAB程序牛顿切线法求方程的近似根（精度为tol）需要自行编制程序.牛顿切线法的MATLAB主程序输入的量：初始值、近似根的误差限tol，近似根的函数值的误差限ftol，迭代次数的最大值gxmax、函数fnq(x)及其导数dfnq(x).输出的量：迭代序列、迭代k次得到的迭代值 (迭代次数超过gxmax时，运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha =和它的偏差的相对误差xdpiancha=的值.根据式（2.12）和已知条件，现提供名为newtonqx.m的M文件：function k,xk,yk,pian

38、cha,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 xdpianchaif (abs(yk)<ftol)&(piancha<tol)|(xdpiancha< tol) k=i-1; xk=x(i);(

39、i-1) xk yk piancha xdpiancha return;endend if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax。') return;end (i-1),xk,yk,piancha,xdpiancha' 用牛顿切线法求方程在的近似根，要求精度.然后判断此方法分别在方程的根处的收敛速度.解 （1）先求方程的近似根.在MATLAB工作窗口输入程序>> k,xk,yk,piancha,xdpiancha=newtonqx(-0.4,0.001, 0.001,100)运行后输出初始值的迭代结果如表 2-10所示

40、.表 2-10kxkyk1-0.516 7-0.076 70.116 70.225 82-0.500 4-0.001 60.016 30.032 63-0.500 0-0.000 00.000 40.000 7迭代次数k = 3,根的近似值xk =-0.500,它的函数值yk =-1.759e-013,偏差piancha =1.712e-007和相对误差xdpiancha =3.423e-007.如果将步骤的迭代结果为k =7；piancha =7.290e-004；xdpiancha =7.295e-004；xk = 0.999；yk =1.590e-006. （2）再讨论收敛速度比较初始值

41、分别是的两组结果，在处迭代3次就得到单根的近似值xk =-0.500；而在处迭代7次才得到二重根的近似值xk单根处的迭代速度比二重根处的迭代速度快很多.这正如前面讨论的结果一样，即若是的单根，则牛顿切线法是平方收敛的；若是的二重根，则牛顿切线法是线性收敛的；我们也可以用阶的定义判断它们的收敛速度.留给读者证明.求的方法及其MATLAB程序用的次根的迭代公式求的近似值（精度为tol）需要自行编制程序.求的次根（当n是偶数时，）的MATLAB主程序输入的量：初始值、根指数、被开方数、的误差限tol，迭代次数的最大值gxmax.输出的量：迭代序列、迭代k次得到的迭代值 (迭代次数达到最大值gxmax

42、时，运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha =|和它的偏差的相对误差xdpiancha=的值.根据求的次根的迭代公式(2.25)和已知条件，现提供名为kai2fang.m的M文件：function k,xk,yk,piancha,xdpiancha,P=kainfang(x0,c,n,tol, gxmax)x(1)=x0;for i=1: gxmaxu(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=pia

43、ncha/( abs(x(i+1)+eps);i=i+1; xk=x(i);yk=fnq(x(i);(i-1),xk,yk,piancha,xdpiancha if (piancha<tol)|(xdpiancha< tol) k=i-1;xk=x(i); return;endendif i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax.')(i-1),xk,yk,piancha,xdpiancha return;endP=(i-1),xk,yk,piancha,xdpiancha' 求，要求精度为.解 本题介绍四种解法.方法1

44、用求的次根（当n是偶数时，）的MATLAB程序计算.取初始值，根指数，被开方数，近似根的精度tol=，迭代的最大次数gxmax=100.在工作区间输入程序>> k,xk,yk,piancha,xdpiancha=kainfang(10,113,2,1e-5,100)运行后输出结果可见,10.630 15，满足精度.方法2 用牛顿迭代公式（2.12）计算.设，则，记，.由牛顿迭代公式得，即取初始值，计算结果列入表 2-12. 表 2-12迭代次数偏差根的近似值10.650 00010.650 00020.019 83610.630 16430.000 01910.630 14640.

45、000 00010.630 146因为，迭代次数=4时，偏差，满足精度，所以，10.630 15.方法3 用牛顿切线法的MATLAB主程序计算.分别建立名为fnq.m和dfnq.m的M文件function y=fnq(x) y=x2-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次，得到精度为的结果 10.630 15.表 2-13kpianchaxdpianchaxkyk1

46、0.650 0000.061 03310.650 0000.422 50020.019 8360.001 86610.630 1640.000 39330.000 0190.000 00210.630 1460.000 00040.000 0000.000 00010.630 1460.000 000方法4 在MATLAB工作空间输入程序>> 113(1/2) 运行后输出经过四舍五入后，得到精度为的结果10.630 15.牛顿切线法的加速及其两种MATLAB程序当是方程的重根时，由牛顿切线公式（2.12）产生的迭代序列收敛到的速度较慢，是线性收敛的.我们要提高收敛速度，可以有选择地

47、采用已知根的重数和未知根的重数两种求重根的修正牛顿切线公式.（一） 已知方程根的重数已知方程根的重数,求重根的修正牛顿切线法的MATLAB主程序输入的量：初始值，方程根的重数，近似根的误差限tol，近似根的函数值误差限ftol，迭代次数的最大值gxmax、函数fnq(x)及其导数dfnq(x)；输出的量：迭代序列、迭代k次得到的迭代值(迭代次数达到最大值gxmax时运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha =|和它的偏差的相对误差xdpiancha=的值. 根据式（2.26）和已知条件，现提供名为newtonxz.m的M文

48、件：function k,piancha,xdpiancha,xk,yk=newtonxz(m,x0,tol,ftol,gxmax)x(1)=x0;for i=1: gxmaxx(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<tol)|(xdpiancha< tol)&(abs(yk)<ftol) k=i-1; xk=x(i);(i-1)

49、piancha xdpiancha xk yk return;endend if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax.') return;endP=(i-1),piancha,xdpiancha,xk,yk' 判断是方程e的几重根？在区间上，分别用牛顿切线法和求重根的修正牛顿切线公式求此根的近似值，使精确到.解 经判断知，是方程的三重根.在MATLAB工作窗口输入程序>> k,piancha,xdpiancha,xk,yk= newtonqx (0.5,1e-4, 1e-4,100)运行后整理结果列入表 2-20.表

50、 2-20kpianchaxdpianchaxkyk10.144 100.404 890.355 900.541 9820.107 410.432 250.248 490.168 2030.077 450.452 800.171 040.051 4640.054 500.467 600.116 550.015 5850.037 690.477 980.078 850.004 6960.025 760.485 140.053 100.001 4070.017 460.490 010.035 630.000 4280.011 770.493 300.023 860.000 1290.007 910

51、.495 520.015 960.000 04100.005 300.497 000.010 660.000 01110.003 540.498 000.007 120.000 00120.002 370.498 660.004 750.000 00130.001 580.499 110.003 170.000 00140.001 050.499 400.002 110.000 00150.000 700.499 580.001 410.000 00160.000 470.499 690.000 940.000 00170.000 310.499 730.000 630.000 00180.0

52、00 210.499 670.000 420.000 00190.000 140.499 440.000 280.000 00200.000 090.498 860.000 190.000 00迭代次数k=20，精确到的根的近似值是xk =0.000 2，其函数值是yk =5.755 8e-011，是一阶（线性）收敛速度.根据求重根的修正牛顿切线公式，在MATLAB工作窗口输入程序 >> k,piancha,xdpiancha,xk,yk= newtonxz(3,0.5,1e-4, 1e-4,100)运行后整理结果得表 2-21.表 2-21kpianchaxdpianchaxky

53、k10.432 306.385 790.067 700.002 9420.066 5457.310 770.001 160.000 0030.001 162673.312 010.000 00-0.000 0040.000 130.996 700.000 130.000 0050.000 13151.535 850.000 00-0.000 0060.000 100.991 440.000 100.000 0070.000 1088.974 840.000 00-8.881 78e-016迭代次数k=7，精确到的根的近似值是xk =1.121 1e-006，piancha = 9.975 3e-005，xdpiancha=88.974 8，其函数值是yk = -8.881 8e-016，是二阶收敛（平方收敛）速度.可见，求重根的修正牛顿切线公式比牛顿切线法收敛速度快得多.（二） 未知方程根的重数未知方程根的重数为，求重根的修正牛顿切线法的MATLAB主程序输入的量：初始值，近似根的误差限to