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

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 第二章 非线性方程（组）的数值解法 2.1 方程（组）的根及其MATLAB命令

2.1.2 求解方程（组）的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)

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

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

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

2.2.1 作图法及其MATLAB程序

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

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

grid, gtext('y=f(x)')

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

⑵ 区间[a,b] 的两个端点的距离 b-a 和步长h的绝对值越小，图形越精确. 作函数y?f(x)在区间 [a,b]上的图形的MATLAB程序二

将y?f(x)化为h(x)?g(x)，其中h(x)和g(x)是两个相等的简单函数.

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

grid,gtext(' y1=h(x),y2=g(x)')

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

2.2.2 逐步搜索法及其MATLAB程序

逐步搜索法的MATLAB主程序

function [k,r]=zhubuss(a,b,h,tol)

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

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

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

7．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 % --- 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程序分别求方程2xsin(cos33?2x2?3x?3?0和

2x)?0在区间[?2,2]上的根的近似值，要求精度是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个是方程2x?2x?3x?3?0的近似根，即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)) 代替，可得到方程

3sin(cos2x)?0在区间[?2,2]上的根的近似值如下

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

1.3310 1.5780 1.7650 1.9200

32

2.3 二分法及其MATLAB程序

2.3.2 二分法的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

8．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 a=x; ya=yx; end

if b-a< abtol , return, end end

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

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

3

解 （1）先用两种方法确定方程x-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的图像.从图像可以看出，此曲线有两个驻点?33都在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)满足f(?2)?f(?1)?0，所以此方程在(-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 0 1 2 3 4 5 6 7 8 9 -2.000 0 -2.000 0 -2.000 0 -1.875 0 -1.812 5 -1.812 5 -1.796 9 -1.796 9 -1.796 9 -1.796 9 -1.000 0 -1.500 0 -1.750 0 -1.750 0 -1.750 0 -1.781 3 -1.781 3 -1.789 1 -1.793 0 -1.794 9 中点xk -1.500 0 -1.750 0 -1.875 0 -1.812 5 -1.781 3 -1.796 9 -1.789 1 -1.793 0 -1.794 9 -1.795 9 bk?ak2 函数值f(ak) -2.000 0 -2.000 0 -2.000 0 -0.716 8 -0.141 8 -0.141 8 -0.004 8 -0.004 8 -0.004 8 -0.004 8 函数值f(bk) 4.000 0 2.125 0 0.390 6 0.390 6 0.390 6 0.129 6 0.129 6 0.062 7 0.029 0 0.012 1 函数值f(xk) 2.125 0 0.390 6 -0.716 8 -0.141 8 0.129 6 -0.004 8 0.062 7 0.029 0 0.012 1 0.003 7 0.500 0 0.250 0 0.125 0 0.062 5 0.031 3 0.015 6 0.007 8 0.003 9 0.002 0 0.001 0 2.4 迭代法及其MATLAB程序

2.4.2 迭代法的MATLAB程序1

迭代法的MATLAB主程序1

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

9．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 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 求方程f(x)?x2?2x?10的一个正根.

解 在MATLAB工作窗口输入程序

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

2运行后输出用迭代公式xk?1?(10?xk)/2的结果

[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相差越来越大，即迭代序列{xk}发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5.

用迭代公式xk?1?10/(xk?2)运行后输出的结果 [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与根x=2.316 624 790 355 40接近，则迭代序列{xk}收敛，但收敛速度较慢，此迭代法较为成功.

用迭代公式xk?1?xk?(xk2*?2xk?10)/(2xk?2)运行后输出的结果

10．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 [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与根x=2.316 624 790 355 40相差无几，则迭代序列{xk}收敛，且收敛速度很快，此迭代法成功.

*

2.4.5

迭代法的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

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

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

xk?1?(xk?1)/3 (k?0,1,2,?).

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.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程序

11．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序

2.5.2 加权迭代法的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

end end

if i>ddmax

disp('迭代次数超过给定的最大值ddmax') k=i-1; xk=x1(i); [(i-1) xk yk]; return; end

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

例2.5.1 求2e?x?x?0在0.85附近的一个近似根，要求精度??10. 解 在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.

?6

2.5.4 艾特肯(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

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

12．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 return; end

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

例2.5.3 用艾特肯加速方法求2e?x?x?0在0.85附近的一个近似根，要求精度

?6??10.

解 在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程序

牛顿切线法的收敛性及其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 用牛顿切线法的局部收敛性判别方程 exsinx?4的近似根时，由下列初始值x0产生的迭代序列是否收敛？

⑴x0?-1; ⑵x0?0; ⑶x0?1; ⑷x0?2; ⑸ x0?5.5；⑹x0?8. 解 在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如下

13．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 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)

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 用牛顿切线法求方程2x?3x?1?0在x0??0.4和0.9的近似根，要求精度??10.

解

在MATLAB工作窗口输入程序

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

14．

?332第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序

2.6.4 求nc的方法及其MATLAB程序

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

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

例2.6.5 求113，要求精度为10解 本题介绍四种解法.

?6.

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

取初始值x0?10，根指数n?2，被开方数c?113，近似根的精度tol=10?5，迭代的最大次数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 可见,113?10.630 15，满足精度10.

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

设

x?113，则x2?5?113'?0，记f(x)?x2?113，f(x)?2x.由牛顿迭代公式得，

xk?1?xk?xk?1132xk2，即xk?1?（xk?21113xk）(其中k?0,1，2，3，?)

取初始值x0?10，计算结果列入表 2-12.

表 2-12 因为，迭代次数k=4时，偏差偏差xk?1?xk 根的近似值xk 迭代次数k ?5?8x4?x3?10，满足精度10，所以，

1 0.650 000 10.650 000 2 0.019 836 10.630 164 113?10.630 15.

3 0.000 019 10.630 146 方法3 用牛顿切线法的MATLAB主程序计

4 0.000 000 10.630 146 算.

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

function y=fnq(x)

15．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 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?5的结果113? 10.630 15.

表 2-13

k 1 2 3 4

piancha 0.650 000 0.019 836 0.000 019 0.000 000 xdpiancha 0.061 033 0.001 866 0.000 002 0.000 000 xk 10.650 000 10.630 164 10.630 146 10.630 146 yk 0.422 500 0.000 393 0.000 000 0.000 000 方法4 在MATLAB工作空间输入程序

>> 113^(1/2)

运行后输出

ans = 10.63014581273465 经过四舍五入后，得到精度为10?5的结果113?10.630 15.

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

（一） 已知方程根的重数m

已知方程f(x)?0根的重数m,求重根的修正牛顿切线法的MATLAB主程序 现提供名为newtonxz.m的M文件：

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)); [(i-1) piancha xdpiancha xk yk]; if ((piancha

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

return; end end

if i>gxmax

disp('请注意：迭代次数超过给定的最大值gxmax.') k=i-1; xk=x(i); yk=fnq(x(i));

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

3x2例2.6.7 判断x?0是方程f(x)?2e?9x?6x?2?0的几重根？在区间

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

??10?4.

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

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

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

表 2-20 k piancha xdpiancha xk yk 16．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

0.144 10 0.107 41 0.077 45 0.054 50 0.037 69 0.025 76 0.017 46 0.011 77 0.007 91 0.005 30 0.003 54 0.002 37 0.001 58 0.001 05 0.000 70 0.000 47 0.000 31 0.000 21 0.000 14 0.000 09 ?40.404 89 0.432 25 0.452 80 0.467 60 0.477 98 0.485 14 0.490 01 0.493 30 0.495 52 0.497 00 0.498 00 0.498 67 0.499 11 0.499 41 0.499 60 0.499 74 0.499 82 0.499 88 0.499 92 0.499 95 0.355 90 0.248 49 0.171 04 0.116 55 0.078 85 0.053 10 0.035 63 0.023 86 0.015 96 0.010 66 0.007 12 0.004 75 0.003 17 0.002 11 0.001 41 0.000 94 0.000 63 0.000 42 0.000 28 0.000 19 0.541 98 0.168 20 0.051 46 0.015 58 0.004 69 0.001 40 0.000 42 0.000 12 0.000 04 0.000 01 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 迭代次数k=20，精确到??10的根的近似值是xk =0.000 2，其函数值是 yk =5.752 0e-011，收敛速度较慢.

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

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

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

表 2-21 k piancha xdpiancha 6.385 79 57.310 85 1.000 78 1.000 00 356.869 17 ?4xk 0.067 70 0.001 16 0.000 00 -0.000 43 0.000 00 0.011 15 0.000 03 0.000 00 yk 0.002 94 0.000 00 0.000 00 -0.000 00 -0.000 00 0.000 01 0.000 00 0.000 00 1 0.432 30 2 0.066 54 4 0.000 43 6 0.011 15 7 0.011 12 3 0.001 16 3443.447 27 5 0.000 43 9228.792 13 8 0.000 03 1638.927 03 迭代次数k=8，精确到??10的根的近似值是xk =1.900 3e-008，其函数值是yk = 4.440 9e-016，piancha = 3.114 5e-005，xdpiancha=1638.927 0，是二阶收敛（平方收敛）.可见，求重根的修正牛顿切线公式比牛顿切线法收敛速度快得多.

（二） 未知方程根的重数

未知方程f(x)?0根的重数为m，求重根的修正牛顿切线法的MATLAB主程序 现提供名为newtonxz1.m的M文件

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);

17．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 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.')

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

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

解 在MATLAB工作窗口输入程序

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

表 2-22

?4k 1 2 3 4

piancha xdpiancha 0.599 23 6.038 57 0.096 98 43.054 85 0.002 25 1 778.327 84 0 0 xk yk -0.099 23 -0.008 18 -0.002 25 -0.000 00 -0.000 00 0 -0.000 00 0 2.7 割线法及其MATLAB程序

2.7.2 割线法的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)

[(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

2.8 抛物线法及其MATLAB程序

18．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 2.8.2 抛物线法的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)

end end

if i>gxmax

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

x1,x2,x3') return; end

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

2x2例2.8.1 用抛物线法求方程f(x)?e?3x?0的一个实根的近似值xk，使精确到

??10.

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

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

1e-4,100) 运行后输出结果

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

19．

?4第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 ans =

Columns 1 through 4

1.00000000000000

-0.39066350562239

Column 5

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

Columns 1 through 4

2.00000000000000

-0.39064638365394

Column 5

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

Columns 1 through 4

3.00000000000000

-0.39064638082059

Column 5

0.00000000000000 k = 3

piancha =

1.712196845282676e-005 xdpiancha =

4.382983989938760e-005 xk =

-0.39064638082059 yk =

2.220446049250313e-016

0.10000000000000 0.20000000000000

0.00933649437761 0.02389906977038

0.00001712196845 0.00004382983990

即，迭代k =3次，得到精度为10?4的根的近似值xk =-0.390 6，其函数值为yk =2.220 4e-016，xk的相邻最小偏差为piancha =1.712e-005，其相对误差为xdpiancha =4.383 0e-005.

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

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

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

function[ci,D,danfan,xddf,hanfan,Xk,Yk]=newtonzu2(X,tol,fto

l,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.y0')

return;

20．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 end end

22??x?y?16,例2.9.1 用迭代公式（2.46）求解方程组 ?

22??x?y?2.解 这里介绍两种求解已知方程组的方法.

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

（1）设f1(x)?x?y?16，

f2(x)?x?y?2，

2222则

?f1(x,y)?x?2x，

?f1(x,y)?y?2y，

?f2(x,y)?x?2x，

?f2(x,y)?y??2y.

（2）作函数f1(x)和f2(x)的图形.输入程序

>>syms x y

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

运行后屏幕显示图形.取两个初始值x0?2,y0?2.

（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文件； （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;ngmax =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]’;

21．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 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?(2,2)T，迭代次数n=6，精度ε=10-6，

x6?3.000000,y6?2.645751与精确解的误差已不超过ε.

2.9.2 求解n元非线性方程组的牛顿法及其MATLAB程序 解n元非线性方程组的牛顿法的MATLAB主程序

现提供名为newtonzu2.m的M文件

function [ci,D,danfan, xddf, hanfan , Xk, Yk]= newtonzun (X,

tol,ftol,gxmax)

Y=Z(X);

for i=1:gxmax

dY=JZ(X);D=det(dY);Xk=X-(dY\\Y')'; 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-(dY\\Y')';Yk=Y;[ci,D,danfan,xddf,hanfan,X,Y]; else

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

if (hanfan

if i>gxmax

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

22??x?y?16,?6例2.9.2 求解方程组 ?2要求精度??10. 2x?y?2,??解 这里介绍两种求解已知方程组的方法.

方法1 用解n元非线性方程组的牛顿法的MATLAB主程序求解 在MATLAB工作窗口输入程序

>> X=[2,2];

[k,D,danfan, xddf, hanfan , Xk, Yk]= newtonzun (X,1e-6,1e-6,100)

运行后输出结果

k=5，D=-63.49803146638493，danfan=3.932201289951168e-011 xddf=9.830503224877920e-012，hanfan=3.336593498083849e-010 Xk = 3.00000000000000 2.64575131106459 Yk = 1.0e-015 *

0 -0.88817841970013

方法2 用解矩阵方程的方法求解 在MATLAB工作窗口输入程序

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

22．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 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('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值’); end end

[k’,x’,y’]

运行后输出结果

ans =

1.00000000000000 2.00000000000000 2.00000000000000 2.00000000000000 3.25000000000000 2.75000000000000 3.00000000000000 3.00961538461538 2.64772727272727 4.00000000000000 3.00001536003932 2.64575204838080 5.00000000000000 3.00000000003932 2.64575131106469 6.00000000000000 3.00000000000000 2.64575131106459

2.9.3 求解n元非线性方程组的拟牛顿法及其MATLAB程序

解n元非线性方程组的拟牛顿法的MATLAB程序

function [k,hanfan,danfan,xddf,Xk,Yk]=ninewton(X01,X02,tol,

ftol,gxmax)

X=[X01;X02];n=length(X); for j=1:gxmax

Y=NN(X); A=NF(X);D=det(A); Xk=X-((A+ones(n)*eps)\\Y’)’; hanfan=norm(Y(2));danfan=norm(Xk(2)-X(2)); xddf=danfan/(norm(Xk(2))+eps);j=j+1;

X=Xk;Y=Z(X);Yk=Y;k=j;warning off MATLAB:singularMatrix if D~=0

Xk=X-(A\\Y’)’;Yk=Y;j=j+1; [j-1,D,danfan, xddf, hanfan ,Xk,

Yk]’

else

disp(‘请注意！迭代方程组（2.59）中的矩阵（2.58）行列式的值等于

零.’)

end

if (hanfan

[j-1,D,danfan, xddf, hanfan ,Xk, Yk]’; warning off MATLAB:singularMatrix return; end end

if(j> gxmax)

error(‘请注意：迭代次数超过给定的最大值gxmax.’); end

23．

本文来源：https://www.bwwdw.com/article/zdj5.html