数值计算课程设计报告材料非线性方程求根

数值计算方法课程设计报告
课程设计名称：数值计算方法课程设计题目：非线性方程求根年级专业：
组员姓名学号：
指导教师：
完成时间：
非线性方程求根
一、问题提出
随着科学技术，生产力经济的发展，在科学与工程计算中存在着大量方程求根问题，例如贷款购房问题，工厂的最佳订货问题等都需要求解一类非线性方程的根，首先根据实际问题列出数学模型，确定变量，给出各个条件及相关函数；然后对建立的模型进行具体分析和研究，选择合适的求解方法；编写函数的程序，用计算机求出方程的解，通过所求解分析具体情况．
求解非线性方程的问题有以下几种基本方法。

二分法简单易行，但收敛较慢，仅有线性收敛速度。

而且该方法不能用于求偶数重根或复根，但可以用来确定迭代法的初始值。

牛顿法是方程求根中常用的一种迭代方法，它除了具有简单迭代法的优点外，还具有二阶收敛速度（在单根邻近处）的特点，但牛顿法对初始值选取比较苛刻（必须充分靠近方程的根），否则牛顿法可能不收敛。

弦截法是牛顿法的一种修改，虽然比牛顿法收敛慢，但因它不需计算函数的导数，故有时宁可用弦截法而不用牛顿法，弦截法也要求初始值必须选取得充分靠近方程的根，否则也可能不收敛。

二、背景分析
代数方程的求根问题是一个古老的数学问题。

理论上，n次代数方程在复数域内一定有 n个根(考虑重数)。

早在16世纪就找到了三次、四次方程的求根公式，但直到19世纪才证明大于等于5次的一般代数方程式不能用代数公式求解，而对于超越方程就复杂的多，如果有解，其解可能是一个或几个，也可能是无穷多个。

一般也不存在根的解析表达式。

因此需要研究数值方法求得满足一定精度要求的根的近似解。

牛顿迭代法是牛顿在17世纪提出的一种求解方程．多数方程不存在求根公式，从而求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要．
而在各种科学和工程计算中往往要用到非线性方程组的求解，而牛顿法又是最基础的迭代法，在各种计算力学、控制工程等领域中发挥了不可代替的作用．而在数值计算中，非线性方程组的求解同样具有重要意义．随着计算机技术的成熟和高速发展，对于非线性方程求根问题出现了大量的数学软件（如MATLAB,SAS,SPSSD等），计算机已经成为工程师应用数学解决工程问题的主要运算工具．同时，工程专业的学生对数学教育的需求重点正在从手工演绎和运算能力的培养转变到结合计算机软件进行建模、求解和论证能力的培养．我们采用Matlab数学软件平台，通过实例比较了二分法、牛顿迭代法、
弦截法三种基本方法的优缺点。

三、 基本算法思想与实现 （1）二分法
单变量函数方程：
f （x ）=0
其中，f(x)在闭区间[a ，b]上连续、单调，且f(a)*f(b)<0,则有函数的介值定理可知，方程f （x ）=0在（a ，b ）区间内有且只有一个解*x ，二分法是通过函数在区间端点的符号来确定*x 所在区域，将有根区间缩小到充分小，从而可以求出满足给定精度的根*x 的近似值。

下面研究二分法的几何意义： 设1a =1,1b =b, 区间[]11,b a ，中点1x =
2
1
1b a +及()1x f ，若()1x f =0，则*x =1x ，若 f(1a )*f(1x )<0，令2a =1a ，2b =1x ，则根*x ∈ [2a ,2b ]中，这样就得到长度缩小一半的有根区间[2a ,2b ]，若 f(1b )*f(1x )<0，令2a =1x ，2b =1b ，则根*x ∈ [2a ,2b ]中，这样就得到长度缩小一半的有根区间[2a ,2b ],即f(2a )f(2b )<0,此时2b -2a =
2
1
1a b -,对有根区间[2a ,2b ]重复上述步骤，即分半求中点，判断中电处符号，则可得长度有缩小一半的有根区间[2a ,2b ], 如图所示：
重复上述过程，第n 步就得到根*
x 的近似序列{}n x 及包含*
x 的区间套，如下：
（1）
...
],....[],[],[2211⊃⊃⊃n n b a b a b a
（2）]
,[,0)()(*n n n n b a x b f a f ∈<
（3）
n a -
n
b =
)
(1121---n n b a =…=
1
2
--n a
b (4) ,2n n n b a x +=
且|*x -n x |≤12
--n a
b (n=1,2,3…..) 显然lim n x ，且n x 以等比数列的收敛速度收敛于*x ，因此用二分法求f （x ）=0的实根*x 可以达到任意指定精度。

（2）牛顿迭代法
设方程f(x)=0在其根*
x 的某个领域U(*
x ，δ)内有一阶连续导数，且f ’(*
x ) ≠0。

求f(x)=0的根*
x ，首先要将f(x)=0转化为等价形式()x x ϕ=，并使ϕ (x)满足不动
点迭代的一般理论。

于是我们令ϕ (x)=x+h(x)f(x),可由ϕ ‘(1x )=0来确定h(x)的结构，根据ϕ’(x)=1+h ’(*
x )f(*
x )+h(*
x )f ’(x1)=1+h(*
x )f ’(*
x )=0可得
h(*x )=-1/f ’(*x ) ，由于f ’(x) ≠0，且f ’(x) 连续，因此当h(x)=-1/f ’(x) 时， h ’(x1)=0，即令ϕ (x)=x-f(x)/f ‘(x)，从而有迭代格式
1+k x = )
(')
(k k k x f x f x -
（k=0,1,2，…..） 由于1x ,2x , 3x …….都在U 领域里，从而当B 比较小时，可用f ’(0x )可近似代替f ’(k x ),1+k x = k x -
)
()
(0x f x f k ，此方法称为牛顿迭代法 下面研究牛顿法的几何意义：
设r 是方程f （x ）=0的根，选取0x 作为的r 初始近似值，经过(0x ,f(0x ))做曲线y=f(x)的切线的方程：y=f(0x )+f ’(0x )(x-0x )，求出L 与x 的交点的横坐标1x =0x -f(0x )/f ’(0x ),称1x 为r 的一次近似值经过点（1x ,f(1x )）做切线y=f(x)的切线，并求出该切线与x 轴的交点横坐标：2x =1x -f(1x )/f ’(1x ),2x 称为r 的二次近似值，重复以上操作可以得到r 的近似值序列。

下述三个定理分别讨论了牛顿法的收敛性质：
定理1：对于方程f(x)=0，设f （x ）在[a ，b]上有二阶连续导数且满足下述条件：
（1）f(a)f(b)<0;
（2）f ’(x)≠0, )(x f ''≠0,对任意的x ∈[a,b];
（3）存在
x ∈[a,b]，使f （
x ）
)
(0x f ''>0，
则由牛顿法产生的迭代序列{}n x 收敛于f(x)=0的根*x ，且
)(2)
()(**2
**1lim x f x f x x x x k k k '''=--+∞
→
定理2：对于方程f(x)=0，设f （x ）在[a ，b]上有二阶连续导数且满足下述条件：
（1）f(a)f(b)<0;
（2）对任意的x ∈[a,b]， f ’(x)≠0, )(x f ''≠0
（3）
)()
(a f a f '<b-a, )()(b f b f '<b-a
则对于任何0x ∈[a,b]，由牛顿法产生的迭代序列{}n x 收敛于f(x)=0的根*
x
定理 3：设*x 是方程f(x)=0的根，在*x 的某个开区间内)(x f ''连续且f ’(x)≠0，则
存在δ>0,当
0x ∈
【*
x
δ-,*x +δ】时，由牛顿迭代法1+k x =
)
(')
(k k k x f x f x -
（k=0,1,2，…..）式产生的序列{}n x 是以不低于二阶的收敛速度收敛到*x .
（3）弦截法
设
k x ,
1
-k x 为方程f(x)=0的两个近似根。

用差商得：f(
k
x )-f(
1
-k x )/
k x -
1
-k x ，
代替牛顿迭代公式中的导数 f ’(
k
x )，于是得到如下的迭代公式：
1+k x =k x -)
()
()()
(11----k k k k k x x x f x f x f 。

下面研究割弦法的几何意义：
经过点（
k
x ,f(
k
x )）及点（
1
-k x ,f(
1
-k x )）两点作割线，其点斜式方程为：
Y=f （k x ）-)()()(11k k k k k x x x x x f x f -----，其零点为X=k x - )
()()()
(11----k k k k k x x x f x f x f 把
X 用
1
+k x 表示即得到迭代格式，它又称为双点弦割法，需要两个初值
此割线与 X 轴交点的横坐标就是新的近似值1
-k x ，所以弦截法又称为割线法，如图所
示。

下面三个定理为弦割法收敛定理：
定理1：设f （x ）在其零点*x 的邻域U （*x ，δ）= [*x -δ，*
x +δ] ( δ>0)内有
二阶连续导数，
0)(*≠'x f ，则当0x ∈U （*x ，δ）时，由割弦法式产生的序列{}n x 收敛于*
x ，且收敛的阶为1.618。

定理2：设)(x f ''在区间[a,b] 上连续，且满足下述三点
（1）f(a)f(b)<0;
（2）对任意的x ∈[a,b]，有f ’(x)≠0, )(x f ''≠0
（3）
)()
(a f a f '≤b-a, )()(b f b f '≤b-a
则对于任意初始
x ，1x ∈[a,b]，由弦割法产生的迭代序列
{}n x 收敛于f(x)=0唯一的根
*x
定理 3：设在其零点*x 的邻域U(*x ,δ)=[*x -δ，*
x +δ](δ>0)内有二阶连续导数，
f ’(x)≠0则当0x ∈ U(*x ,δ)时，由弦割1+k x =k x -)
()
()()
(11----k k k k k x x x f x f x f
式产生的序列
{}n x 收敛于*x ，且收敛的阶为1.618。

四、 具体应用实例分析
求解033)(2
3
=--+=x x x x f 在5.1附近的根。

（1）二分法 建立erfen -M 文件：
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.'), retur n 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;
else
if yb*yx>0
b=x;
yb=yx;
else
a=x;
ya=yx;
end
if b-a<abtol , retur n,
end
end
k=max1; x; wuca; yx=fun(x);
建立FUN函数文件：functiony=fun(x) y=x.^3+x.^2-3*x-3; 画图：
>> x=[-10:0.1:10]; >> y=fun(x);
>>plot(x,y);
由图，我们选取区间[-6，6]
输入程序：
[k,x,wuca,yx]=erfen(-6,6,0.001)
运行结果：
k =13；
x =1.7322；
wuca =7.3242e-004；
yx =0.0012
（2）牛顿迭代法
建立newtonqx -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) xkykpianchaxdpiancha]
if (abs(yk)<ftol)&((piancha<tol)|(xdpiancha<tol))
k=i-1;
xk=x(i);[(i-1) xkykpianchaxdpiancha]
return;
end
end
if i>gxmax
disp('请注意：迭代次数超过给定的最大值gxmax。

')
k=i-1;
xk=x(i);
[(i-1) xkykpianchaxdpiancha]
return;
end
[(i-1),xk,yk,piancha,xdpiancha]';
建立FNQ原函数文件：
function y=fnq(x)
y=x.^3+x.^2-3*x-3;
建立DFNQ导函数文件：
function y=dfnq(x)
y=3*x.^2+2*x-3;
输入程序：[k,xk,yk,piancha,xdpiancha]=newtonqx(1.5,0.001,0.001,20)
输出结果：
k =3；
xk =1.7321；
yk =1.0617e-005；
piancha =0.0013；
xdpiancha =7.5560e-004
（3）弦截法
建立gexian文件：
function[k,piancha,xdpiancha,xk,yk]=gexian(x01,x02,tol,ftol,gxmax) x(1)=x01;
x(2)=x02;
fori=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) pianchaxdpianchaxkyk]
if (abs(yk)<ftol)&(( piancha<tol)|(xdpiancha<tol))
k=i-2;
xk=x(i);
yk=fnq(x(i));
[(i-2) pianchaxdpianchaxkyk]；
Return;
end
end
if i>gxmax
disp('请注意：迭代次数超过给定的最大值gxmax.')
k=i-2;
xk=x(i);
yk=fnq(x(i));
return;
end
建立FNQ函数文件：
function y=fnq(x)
y=x.^3+x.^2-3*x-3;
输入程序：（取1.5附近的初始值1.4，1.6两个）
[k,piancha,xdpiancha,xk,yk]=gexian (1.4,1.6,0.001,0.001,20)
输出结果：
k =4；
piancha =5.1561e-004；
xdpiancha =2.9769e-004；
xk =1.7321
yk =3.6710e-005
五、设计总结
根据二分法求解非线性方程根的原理，将所求方程根所在的区间平分为两个小区间，在判断根属于哪个小区间；把有根的小区间再平分为二，再判断根所在的更小的区间，对分；重复这一过程，最后求出所要的近似值。

当所分的小区间的间距越小的时候，得出的方程根结果就越精确，其原因就是所分的小区间间距越小，则就越接近方程等于0的根。

所以最后的结果的精度越高，得到的误差越小；而对于简单迭代法，只有在满足一定条件的情况下，才能求解出在区间上有唯一根，使迭代序列收敛于。

根据牛顿迭代法的原理，求解出非线性方程根的结果可以看出，牛顿迭代法具有平方收敛的速度，所以在迭代过程中只要迭代几次就会得到比较精确的解，并不像简单迭代法，需要迭代多次才能解出较为精确的结果，但是用牛顿迭代法求解时选定的初值要接近方程的解，否则可能得不到收敛的结果。

同时，牛顿迭代法计算量也会相对较大些。

单点弦截法，用选定的两个初值点所对应的函数值连接作弦，用此弦与轴的交点横坐标作为方程根的近似值。

按此方法进行迭代计算，直到满足精度要求为止
六、参考文献
[1]奚梅成。

数值分析方法【M】。

合肥：中国科技技术大学出版社，2007.
[2]薛毅。

数值分析与实验【M】。

北京：北京理工大学出版社，2005.
[3]汪卉琴，刘目楼。

数值分析【M】。

北京：冶金工业出版社，2004.
[4]丁丽娟，程杞元。

数值计算方法【M】。

北京：北京理工大学出版社，2005.
[5]薛定宇，陈阳泉。

高等应用数学问题的MATLAB求解。

北京：清华大学出版社，2008.
七、心得体会
（1）韩建：这次数值分析课程设计我们受益匪浅。

在刚开始拿到题目时，我们提出很多问题，4个人各有自己的想法，产生了分歧。

但经过这一周的时间和体验下来，我们学到的不仅是课本知识，还有团队和合作精神。

现在想来，也许学校安排的课程设计有着它更深层的意义，它不仅仅让我们综合那些理论知识来运用到设计和创新，还让我们知道
了一个团队凝聚在一起是所能发挥出的巨大潜能！
在这次课程设计中，我们运用到了以前所学的专业课知识，牛顿迭代法、matlab汇编语言等。

虽然过去都将这些知识用于解题中，未有独立应用过它们，但在学习的过程中带着问题去学我发现效率很高。

设计过程，好比我们的成长历程，常有一些不如意，难免会遇到各种各样的问题。

这也激发了我今后努力学习的兴趣，通过这次设计，我懂得了学习的重要性，了解到理论知识与实践结合的重要意义，学会了坚持、耐心和努力，这将为自己今后的学习和工作打下基础。

（2）高育坤：此次实验过程中学到了许多东西，对于MATLAB的运用也比以前熟练了一些，在编程过程中也考虑了一些题目以外的因素，虽然还有不足之处，不过我相信只要不断努力一定会有所进步，通过几次实验我也明白了一个道理，对于编程不要一看到就害怕，一步步的去做，一点点的去实现算法，编程最难的是算法，题目给出了算法我们需要的就是将数学公式转化为编程语言
（3）李婧：通过此次课程设计我们可以知道计算机在现代生活中的应用已经如此普及，尤其是在数学计算当中，Matlab软件更是发挥了不可替代的作用．Matlab以其强大的功能，方便了当今数值计算，数学教程，及工程计算等众多领域．我们在以Matlab软件为平台的基础上，给出了非线性方程的一般解法，非线性方程的求解有二分法，牛顿迭代法等．二分法的优点是算法简单，且总是收敛的，但由于二分法的收敛速度太慢，故一般不单独将其用于求根，只用其为根求得一个较好的近似．本文主要介绍了牛顿迭代法及其在现实生活中的应用．牛顿迭代法为平方收敛，故其收敛速度较快，但对初值的选取需要谨慎，如果初值选取错误，则可能导致方程迭代发散，最终不能求解出正确解．在计算一些对精度要求特别苛刻时，最好给出较高的精度输入及输出，防止因为精度问题导致误差过大，最终影响结果．
（4）王冬妮：二分法的优点是计算简单，方法可靠，误差容易估计，只要求连续，且总是收敛的，因此对函数的性质要求较低。

它的缺点是不能求偶数重根，也不能求复根，且收敛较慢。

故一般不单独将其用于求根，只用其为根求得一个较好的近似值。

牛顿迭代法是多项式求根的一种效率很高的算法，收敛速度快（对单根）。

算法简单是迭代法中较好者，但是它有两个缺点：第一每次只能求出一个ε－根，求其它根时若采用降次处理又会产生精度降低的问题。

第二有时会遇到由于初始点选择不当而使算法失效。

牛顿法和割弦法都是先将f(x)线性化，然后求根，但线性化的方式不同：从分析的角度说，牛顿法是在根 *x邻近点处的切线函数作为f(x)的近似，而割弦法是在*x邻近用f(x)的一次插值函数作为f(x)的近似函数，它们本质的区别在于，牛顿迭代法在计算时，只用到前一步的值，弦截法需要用两个猜测值，，因此使用这种方法必须先给出
两个初始值 , 。

根据以上三种方法的优缺点，我们在使用非线性方程求根时，应根据实际方程选出一种或者多种方法进行综合求解，以便快速，便捷的求出最佳精确值。

八、附录
（1）二分法的程序
function [c,err,yc] = bisect(f1,a,b,delta)
% f1是所要求解的函数。

% a和b分别为有根区间的左右限。

% delta是允许的误差界。

% c为所求近似解。

% yc为函数f在c上的值。

% err是c的误差估计。

ya = feval('f1',a);
yb = feval('f',b);
if yb == 0,
c = b；
return
end
if ya*yb>0,
disp('(a,b)不是有根区间');
return
end
max1 = 1 + round((log(b-a) - log(delta))/log(2));
for k = 1:max1
c = (a + b)/2;
yc = feval('f',c);
if yc == 0
a=c;
b=c;
return
elseif yb*yc>0
b = c;
yb = yc;
else
a = c;
ya = yc;
end
if (b-a) < delta,return,end
end
k;
c = (a + b)/2;
err = abs(b-a);
yc = feval('f',c);
****************************************************************** （2）牛顿法的程序
function [p1,err,k,y] = newton(f1,df1,p0,delta,max1)
% f1是非线性函数。

% df1是f1的微商。

% p0是初始值。

% delta是给定允许误差。

% max1是迭代的最大次数。

% p1是牛顿法求得的方程的近似解。

% err是p0的误差估计。

% k是迭代次数。

% y = f(p1)
p0, feval('f1',p0)
for k = 1:max1
p1 = p0 - feval('f1', p0)/feval('df1', p0);
err = abs(p1-p0);
p0 = p1;
p1, err, k, y = feval('f1', p1)
if (err < delta) | (y == 0),
break,
end
p1, err, k, y = feval('f1', p1)
end
************************************************************** （3）弦截法的程序
function [p1, err, k, y] = secant(f1, p0, p1, delta, max1)
% f1是给定的非线性函数。

% p0,p1为初始值。

% delta为给定误差界。

% max1是迭代次数的上限。

% p1为所求得的方程的近似解。

% err为p1-p0的绝对值。

% k为所需的迭代次数。

% y＝f(p1)
p0, p1, feval('f1', p0), feval('f1',p1), k = 0;
for k = 1:max1
p2 = p1 - feval('f1',p1)*(p1-p0)/(feval('f1',p1) - feval('f1042',p0));
err = abs(p2-p1);
p0 = p1;
p1 = p2;
p1, err, k, y = feval('f1', p1)；
if (err < delta) | (y == 0),
break,
end
end
******************************************************************。