牛顿迭代法收敛及算法实验

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 （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

2.6.3 牛顿切线法的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)

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]';

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

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

解

在MATLAB 工作窗口输入程序

>> [k,xk,yk,piancha,xdpiancha]=newtonqx(-0.4,0.001, 0.001,100) 运行后输出初始值4.00-=x 的迭代结果.

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

c 的方法及其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

k=i-1;xk=x(i); yk=fnq(x(i));

return ;