第二章非线性方程（组）的数值解法.doc

7.1 方程（组）的根及其MATLAB命令 7.1.1 求解方程（组）的solve命令 求方程f(x)=q(x)的根可以用MATLAB命令: >> x=solve('方程f(x)=q(x)','待求符号变量x') 求方程组fi(x1,…,xn)=qi(x1,…,xn) (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) 7.1.2 求解方程（组）的fsolve命令 fsolve的调用格式： X=fsolve(F,X0) 7.2 搜索根的方法及其MATLAB程序 7.2.1 作图法及其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程序二 将化为，其中是两个相等的简单函数. x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2) grid,gtext(' y1=h(x),y2=g(x)') 说明：此程序在MATLAB的工作区输入，运行后即可出现函数的图形.两图形交点的横坐标即为所要求的根的近似值. 7.2.2 逐步搜索法及其MATLAB程序 逐步搜索法的MATLAB主程序 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); 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); end end 例7.2.1 用逐步搜索法的MATLAB程序分别求方程和在区间上的根的近似值，要求精度是0.000 1. 解 将逐步搜索法的MATLAB程序保存名为zhubuss.m的M文件. 建立M文件funs.m function y=funs(x) y=2.*x.^3+2.*x.^2-3.*x-3 在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)) 代替，可得到方程在区间上的根的近似值如下 r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.9230 1.3310 1.5780 1.7650 1.9200 7.3 二分法及其MATLAB程序 7.3.1二分法的MATLAB程序 二分法的MATLAB主程序 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.'), 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); 例7.3.1 确定方程x3-x+4=0的实根的分布情况，并用二分法求在开区间 (-2,-1)内的实根的近似值，要求精度为0.001. 解 （1）先用两种方法确定方程x3-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)=x3-x+4的图像.从图像可以看出，此曲线有两个驻点都在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)内有一个实根. （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-3 次数k 左端点ak 右端点bk 中点xk 函数值f(ak) 函数值f(bk) 函数值f(xk) 0 -2.000 0 -1.000 0 -1.500 0 0.500 0 -2.000 0 4.000 0 2.125 0 1 -2.000 0 -1.500 0 -1.750 0 0.250 0 -2.000 0 2.125 0 0.390 6 2 -2.000 0 -1.750 0 -1.875 0 0.125 0 -2.000 0 0.390 6 -0.716 8 3 -1.875 0 -1.750 0 -1.812 5 0.062 5 -0.716 8 0.390 6 -0.141 8 4 -1.812 5 -1.750 0 -1.781 3 0.031 3 -0.141 8 0.390 6 0.129 6 5 -1.812 5 -1.781 3 -1.796 9 0.015 6 -0.141 8 0.129 6 -0.004 8 6 -1.796 9 -1.781 3 -1.789 1 0.007 8 -0.004 8 0.129 6 0.062 7 7 -1.796 9 -1.789 1 -1.793 0 0.003 9 -0.004 8 0.062 7 0.029 0 8 -1.796 9 -1.793 0 -1.794 9 0.002 0 -0.004 8 0.029 0 0.012 1 9 -1.796 9 -1.794 9 -1.795 9 0.001 0 -0.004 8 0.012 1 0.003 7 7.4 迭代法及其MATLAB程序 7.4.1迭代法的MATLAB程序1 迭代法的MATLAB主程序1 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]'; 例7.4.1 求方程的一个正根. 解 在MATLAB工作窗口输入程序 >> [k,piancha,xdpiancha,xk]= diedai1(2,5) 运行后输出用迭代公式的结果 [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 由以上运行后输出的迭代序列与根=2.316 624 790 355 40相差越来越大，即迭代序列发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5. 用迭代公式运行后输出的结果 [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次的迭代值xk =2.331 460 674 157 30与根=2.316 624 790 355 40接近，则迭代序列收敛，但收敛速度较慢，此迭代法较为成功. 用迭代公式运行后输出的结果 [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次的迭代值xk=2.316 624 790 355 40与根=2.316 624 790 355 40相差无几，则迭代序列收敛，且收敛速度很快，此迭代法成功. 7.4.2迭代法的MATLAB程序2 迭代法的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<tol)|(xdpiancha< tol) k=i-1; xk=x(i); return; end end if i>ddmax disp('迭代次数超过给定的最大值ddmax') k=i-1; xk=x(i);yk=fun(x(i)); [(i-1) piancha xdpiancha xk yk]; return; end P=[(i-1),piancha,xdpiancha,xk,yk]'; 例7.4.2 求x5-3x+1=0在0.3附近的根，精确到4位小数. 解 ⑴ 构造迭代公式   . ⑵利用迭代法的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. 7.5 迭代过程的加速方法及其MATLAB程序 7.5.1加权迭代法的MATLAB程序 加权迭代法的MATLAB主程序 现提供名为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; xk=x1(i); return; end end if i>ddmax disp('迭代次数超过给定的最大值ddmax') k=i-1; xk=x1(i); [(i-1) xk yk]; return; end P=[(i-1),xk,yk]'; 例7.5.1 求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.85260550406236 3.00000000000000 0.85260550207699 0.85260550208255 k =3； xk =0.852606； yk = 0.852606. 7.5.2 艾特肯(Aitken)加速方法的MATLAB程序 艾特肯(Aitken)加速方法的MATLAB主程序 现提供名为Aitken.m的M文件如下： 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<tol)|(xdpiancha<tol) k=i-1; xk=x(i); yk=fun(x(i)); m=[0,1:i-1]; p=[m',x1',x2',x']; return; end end if i>ddmax disp('迭代次数超过给定的最大值ddmax') k=i-1; xk=x(i); yk=fun(x(i)); m=[0,1:i-1]; p=[m',x1',x2',x']; return; end m=[0,1:i-1]; p=[m',x1',x2',x']; 例7.5.3 用艾特肯加速方法求2e在0.85附近的一个近似根，要求精度. 解 在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 7.6 牛顿(Newton)切线法及其MATLAB程序 7.6.1 牛顿切线法的收敛性及其MATLAB程序 牛顿切线法的收敛性及其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]'; 例 7.6.1 用牛顿切线法的局部收敛性判别方程 e的近似根时，由下列初始值产生的迭代序列是否收敛？ ⑴ ⑵ ⑶ ⑷ ⑸ ；⑹. 解 在MATLAB工作窗口输入程序 >> [y,f]=newjushou(-1) 运行后输出结果 请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值 y = 139.5644 f = 4.3096 （2）输入程序 >> [y,f]=newjushou(0) 运行后输出结果 请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)的值f y = 8.0000 f = 4 （3）输入程序 >> [y,f]=newjushou(1) 运行后输出结果 恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下 y = 0.3566 f = 1.7126 （4）输入程序： >> [y,f]=newjushou(2) 运行后输出结果 请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值 y = 1.2593 f = -2.7188 （5）输入程序 >> [y,f]=newjushou(5.5) 运行后输出结果 请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值 y = 1.0447e+005 f = 176.6400 （6）输入程序 >> [y,f]=newjushou(8) 运行后输出结果 恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下 y = 0.4038 f = -2.9452e+003 7.6.2 牛顿切线法的MATLAB程序 牛顿切线法的MATLAB主程序 现提供名为newtonqx.m的M文件： 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)<ftol)&((piancha<tol)|(xdpiancha< tol)) k=i-1; xk=x(i);[(i-1) xk yk piancha xdpiancha] return; end end if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax。') k=i-1; xk=x(i);[(i-1) xk yk piancha xdpiancha] return; end [(i-1),xk,yk,piancha,xdpiancha]'; 例 7.6.2 用牛顿切线法求方程在的近似根，要求精度. 解 在MATLAB工作窗口输入程序 >> [k,xk,yk,piancha,xdpiancha]=newtonqx(-0.4,0.001, 0.001,100) 运行后输出初始值的迭代结果. 7.6.3 求的方法及其MATLAB程序 求的方法及其MATLAB主程序 现提供名为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<tol)|(xdpiancha< tol) k=i-1;xk=x(i); yk=fnq(x(i)); return; end end if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax.') k=i-1;xk=x(i); yk=fnq(x(i)); [(i-1),xk,yk,piancha,xdpiancha] return; end P=[(i-1),xk,yk,piancha,xdpiancha]'; 例7.6.3 求，要求精度为. 解 本题介绍四种解法. 方法1 用求的次根（当n是偶数时，）的MATLAB程序计算. 取初始值，根指数，被开方数，近似根的精度tol=，迭代的最大次数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 可见,10.630 15，满足精度. 方法2 用牛顿迭代公式计算. 设，则，记，.由牛顿迭代公式得，，即 取初始值，计算结果列入表 迭代次数 偏差 根的近似值 1 0.650 000 10.650 000 2 0.019 836 10.630 164 3 0.000 019 10.630 146 4 0.000 000 10.630 146 因为，迭代次数=4时，偏差，满足精度，所以，10.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次，得到精度为的结果 10.630 15. k piancha xdpiancha xk yk 1 0.650 000 0.061 033 10.650 000 0.422 500 2 0.019 836 0.001 866 10.630 164 0.000 393 3 0.000 019 0.000 002 10.630 146 0.000 000 4 0.000 000 0.000 000 10.630 146 0.000 000 方法4 在MATLAB工作空间输入程序 >> 113^(1/2) 运行后输出 ans = 10.63014581273465 经过四舍五入后，得到精度为的结果10.630 15. 7.7 割线法及其MATLAB程序 7.7.1 割线法的MATLAB程序 割线法的MATLAB主程序 现提供名为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)<ftol)&(( piancha <tol)|(xdpiancha< tol)) k=i-2; xk=x(i);yk=fnq(x(i)); [(i-2) piancha xdpiancha xk yk]； return; end end if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax.') k=i-2; xk=x(i);yk=fnq(x(i)); return; end 7.8 抛物线法及其MATLAB程序 7.8.1 抛物线法的MATLAB程序 抛物线法的MATLAB主程序 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)<ftol)&(( piancha <tol)|(xdpiancha< tol)) k=i; X(3)= xk; Y(3)=fnq(X(3