数值分析课程设计
2013/2014第一学期
数值分析课程设计
设计题目：
非线性方程的数值解法及MATLAB 解法 专业
学号 姓名 学号 姓名 学号 姓名
成绩 指导老师
摘 要
本论文分析总结了非线性方程的求解的几种方法，主要介绍非线性方程的数值解法，是直接从方程出发，逐步缩小根的存在区间，或逐步将根的近似值精确化，直到满足问题对精度的要求。分别介绍了二分法，黄金分割法，迭代法，Newton 法，弦截法，抛物法，还详细编写注释了这几种方法的Matlab 函数代码。并通过具体的例子对他们做出分析，通过各种方法的对照比较得出各个算法的优缺点。
关键词：非线性方程 二分法 黄金分割法 迭代法 弦截法 Matlab
问 题 提 出
1.5]1. 分别用二分法，黄金分割法求方程f (x ) =x 3-x -1=0在区间[1.0,
内的一个实根，要求准确到小数点后第2位。（书P214）
2. 用迭代法求方程f (x ) =x 3-x -1=0在x 0=1.5附近的根x *。（书P215） 3. 分别用牛顿法，弦截法，抛物法求方程f (x ) =xe x -1=0。（书P223，P229）
1 二分法的基本原理及MATLAB 实现
二分法的基本思想：
介值定理：若函数在区间[a , b ]上连续，且f (a ) f (x )
设区间[a , b ]是有根区间，令a 0=a , b 0=b , 取区间中点 x 0=
a 0+b 0
2
。
考查中点x 0的函数值情况，若f (a 0) f (x 0)
此时令a 1=a 0, b 1=x 0；否则f (a 0) f (x 0) >0（不考虑等号情况，x *=[a 0, b 0]，
否则f (x 0) =0，已得到方程的根），因此得到f (a 0) 则x *∈[x 0, b 0]，f (x 0)
此时令a 1=x 0, b 1=b 0，再取x 1=
a 1+b 1
2
。
设当前的有根区间[a k , b k ], 取x k =
a k +b k
2
，若f (a k ) f (x k )
a k +1=a k , b k +1=x k ；否则令a k +1=x k , b k +1=b k ；再取x k =
入下一轮计算。 二分法的几何意义：
逐步将区间对分，得到根的近似值。 二分法的收敛性质：
a k +1+b k +1
2
，转
设初始有根区间为[a , b ]，x *是方程的根，x k 为第k 次区间[a k , b k ]的中点（初始区间记为0次）则
k -x *≤
b k -a k
2
=
b k -1-a k -1
2
2
= =
b -a
2
k +1
上式给出了算法的终止条件，当
b k -a k
2
给定的精度要求。上式同时表明。当k →∞有x k -x *→0。因此序列的收敛性与初始区间[a , b ]无关，二分法是大范围收敛的。 下面是二分法的MATLAB 函数代码： %二分法求解非线性方程 function root=HalfInterval(f,a,b,eps) if(nargin==3) eps=1.0e-4; end
f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end
if(f2==0) root=b; end if(f1*f2>0)
disp('两端点函数值乘积大于0！'); return; else
root=FindRoots(f,a,b,eps); end
function r=FindRoots(f,a,b,eps) f_1=subs(sym(f),findsym(sym(f)),a); %f_2=subs(sym(f),findsym(sym(f)),b); mf=subs(sym(f),findsym(sym(f)),(a+b)/2); if(f_1*mf>0) t=(a+b)/2;
r=FindRoots(f,t,b,eps); else
if(f_1*mf==0) r=(a+b)/2; else
if(abs(b-a)
s=(a+b)/2;
r=FindRoots(f,a,s,eps); end end end 解：
在Matlab 命令窗口输入：
root=HalfInterval('x^3-x-1',1.0,1.5,1.0e-2) 计算出结果为： root = 1.3252
2 黄金分割法的基本原理及MATLAB 实现
二分法是把区间的长度减半，黄金分割法是把区间逐步缩短为前一次的0.618倍，其求解步骤如下：
1）
设t 1=a +(1-0. 618) *(b -a ) , t 2=a +(1-0. 618) *(b -a ) ，且f 1=f (t 1) , f 2=f (t 2) .
2）
如果1-t 2
t 1+t 2
2
；否则转到3）。
3）
如果f 1*f 2
f 1*f (a )
下面是黄金分割法的MATLAB 函数代码： %黄金分割法求解非线性方程 function root=GoldenSection(f,a,b,eps) if(nargin==3) eps=1.0e-4; end
f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0)
root=b; end if(f1*f2>0)
disp('两端点函数值乘积大于0！'); return; else
t1=a+(b-a)*0.382; t2=a+(b-a)*0.618;
f_1=subs(sym(f),findsym(sym(f)),t1); f_2=subs(sym(f),findsym(sym(f)),t2); tol=abs(t1-t2);
while(tol>eps) if(f_1*f_2
fa=subs(sym(f),findsym(sym(f)),a); if(f_1*fa>0) a=t2; else b=t1; end end
t1=a+(b-a)*0.382; t2=a+(b-a)*0.618;
f_1=subs(sym(f),findsym(sym(f)),t1); f_2=subs(sym(f),findsym(sym(f)),t2); tol=abs(t2-t1); end
%精度控制
root=(t1+t2)/2; %输出根 end 解：
在Matlab 命令窗口输入：
root= GoldenSection ('x^3-x-1',1.0,1.5,1.0e-2) 计算出结果为： root = 1.3229
3 迭代法的基本原理及MATLAB 实现
迭代法的基本原理
其原理是构造一个迭代公式，反复用它得出一个逐次逼近方程的数列，数列
中的每个元素方程根的近似值，只是精度不同而已。
迭代法求解方程f (x ) =0时，先把方程等价地变换成形式
f (x ) =x -g (x ) =0，移项得出x =g (x ) ，若函数g (x ) 连续，则称式x =g (x ) 为迭代函数。用它可构造出迭代公式
x k +1=g (x k ), k =0, 1, 2,
从初始值x 0出发，便可得出迭代序列
{x k }=x , x , x , , x k ,
1
2
如果迭代序列收敛，且收敛于x *，则由式x k +1=g (x k ), k =0, 1, 2, 有
lim (g (x k ) -x k +1) =(g (x *) -x *) =f (x *) =0
k →∞
可见，x *便是方程式f (x ) =0的根。 迭代法的几何意义
解方程f (x ) =0可以等价地变换成求解x =g (x ) ，这就等于求曲线y =x 和y =g (x ) 交点P *的坐标x *。求迭代序列，就等于从图中点x 0出发，由函数
y =g (x 0) 得出y =P 0，代入函数y =x 中得到Q 1，再把Q 1的x 坐标x 1代入方程
y =g (x ) 得出P 1，如此继续下去，便可在曲线y =g (x ) 上得到一系列的点P 0, P 1, , P k , ，这些点的x 坐标便是迭代数列x 1, x 2, , x k , ，它趋向于方程式f (x ) =0的根x *，数列的元素就是方程根的近似值，数列的收敛等价于曲线
y =x 和y =g (x ) 能相交于一点。
迭代公式收敛定理
方程x =g (x ) 在(a , b )内有根x *，如果 1） 当x ∈[a , b ]时，g (x ) ∈[a , b ]。
2） g (x ) 可导，且存在正数q
1． 方程x =g (x ) 在(a , b )内有唯一的根x *。
2． 迭代公式x k +1=g (x k ) 对(a , b )内的任意初始近似根x 0均收敛于x *。 3． 近视根x k 的误差估计公式为
x -x k ≤
*
q k
1-q
x 1-x 0
由该定理知，将方程f (x ) =0转化为等价形式x =g (x ) 时，选择和构造什么样的迭代函数g (x ) 非常重要，只有当它满足一定条件时，迭代序列才收敛于方程的根x *。
下面是迭代法的MATLAB 函数代码： %迭代法求解非线性方程
function [p0,k,err,p]=Iteration(g,p0,tol,max1) %g是给定的迭代函数 %p0是给定的初始值
%max1是所允许的最大迭代次数 %k所进行的迭代次数加1 %p是不动点的近似值
%err是误差 %p(p1,p2……pn) P(1)=p0; for k=2:max1
P(k)=feval('g',P(k-1)); k,err=abs(P(k)-P(k-1)) p=P(k); if(err
disp('超过了最大的迭代数量'); end end P 解：
先用M-文件写一个名为g.m 的函数。 %迭代法实例函数 function y=g(x) y=x^3-x-1;
在Matlab 命令窗口输入： Iteration('g',1.5,10^(-4),10)
由运行结果可知：该迭代过程是发散的，所以没有意义，应重新选取迭代公式。
4 牛顿法的基本原理及MATLAB 实现
牛顿法的基本原理
将非线性方程线性化，以线性方程的解逐步逼近非线性方程的解。 把函数f (x ) 在某一初始值x 0点附近展开成泰勒级数，有
f (x ) =f (x 0) +(x -x 0) f ' (x 0) +(x -x 0)
2
f , , (x 0)
2!
+
取其线性部分，近似的代替函数f (x ) 可得方程的近似式
f (x ) ≈f (x 0) +(x -x 0) f ' (x 0) =0
(x ) ≠0，解该近似方程可得 设f ’
x 1=x 0-
f (x 0)
f ' (x 0)
x 1可作为方程式的近似解，重复以上过程，得迭代公式
x k +1=x k -
按上式可求方程式的近似解。
f (x k )
f ' (x k )
牛顿法的收敛性
(x ) ≠0，f ’’(x ) 连续且不变号，则只要选取的初在有根区间[a , b ]上，f ’
始值近似根x 0满足f (x 0) f ’' (x ) >0，牛顿法必定收敛。
用牛顿迭代公式在某次算得的误差与上次误差的平方成正比。可见，牛顿迭代公式的收敛速度很快。但计算实践表明，当初始值不够好时，牛顿法可能发散。一般可由问题的实际背景来预测或由对分区间法求的较好的初值。 下面是牛顿法的MATLAB 函数代码： %牛顿法求解非线性方程
function [p1,err,k,y]=Newton(f,df,p0,delta,max1) %f是非线性函数
%df是f 的微商 %p0是初始值 %delta是给定允许误差 %max1是迭代的最大次数 %p1是牛顿法求得的方程的近似值 %err是p0的误差估计 %k是迭代次数 %y=f(p1) p0, feval('f',p0) for k=1:max1
p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; p1, err, k,
y=feval('f',p1) if(err
首先在Matlab 输入
fplot('[x*exp(x)-1]',[-1,1]);grid
回车得到下图，可知函数f (x ) 与x 轴在x =0. 5附近有根，取处值x 0=0. 5。
再用M-文件写一个名为f.m 的函数。 %牛顿法, 弦截法实例函数f function y=f(x) y=x*exp(x)-1; 名为df.m 的函数 %牛顿法实例函数df function y=df(x) y=(1+x)*exp(x);
在Matlab 命令窗口输入： Newton('f','df',0.5,10^(-4),10) 运行得 ans =
0.56714
5 弦截法的基本原理及MATLAB 实现
弦截法的基本原理
给定非线性方程f (x ) =0，选定曲线y =f (x ) 上的两个点P 0(x 0, f (x 0)) ，
P 1(x 1, f (x 1)) ，过这两个点做一条直线P 0P 1，则直线方程
y =f (x 0) +
f (x 1) -f (x 0)
(x -x 1) 。
x 1+x 0
当f (x 0) =f (x 1) 时，直线P 0P 1与x 轴的交点为x 2=x 1+
f (x 1)(x 1-x 0)
，这时
f (x 1) -f (x 0)
用x 2作为曲线y =f (x ) 与x 轴交点的近似值，显然
x 2-x *
这里的x *为f (x ) =0的精确解。然后用P 1(x 1, f (x 1)) ，P 2(x 2, f (x 2)) , 构造直线
)
P 1P 2。重复上述步骤，就可以求出x 3。
如此进行下去，可得到迭代格式
x k +1=x k +
f (x k )(x k -x k -1)
f (x k ) -f (x k -1)
f (x k ) -f (x k -1) f (x k )
取代牛顿公式x k +1=x k -中的
f (x k ) -f (x k -1) f ' (x k )
上式实际上就是用均差
微商f ' (x ) 的结果，所以弦截法可以被看成是牛顿法的一种变形。 下面是弦截法的MATLAB 函数代码： %弦截法求解非线性方程
function [p1,err,k,y]=Secant(f,p0,p1,delta,max1) %f是给定的非线性函数 %p0,p1为初始值 %delta为给定误差界 %max1是迭代次数的上限 %p1为所求得的方程的近似解 %err为p1-p0的绝对值 %k为所需要的迭代次数 %y=f(p1)
k=0,p0,p1,feval('f',p0),feval('f',p1) for k=1:max1
p2=p1-feval('f',p1)*(p1-p0)/(feval('f',p1)-feval('f',p0)); err=abs(p2-p1); p0=p1; p1=p2;
k,p1,err,y=feval('f',p1)
if(err
在Matlab 命令窗口输入： Secant('f',0.5,0.6,10^(-4),10) 回车运行得 ans =
0.56714
6 抛物法的基本原理及MATLAB 实现
抛物法的基本原理
设已知方程f (x ) =0的3个近似根为x k , x k -1, x k -2，我们以由这3个点为节点构造出的二次插值多项式P (x ) 的一个零点x k -1作为新的近似根，这样确定的迭代过程称为抛物法。
在几何图形上，这种方法的基本思想是用抛物线与x 轴的交点x k -1作为所求根的近似位置。
现在推导抛物法的计算公式。 插值多项式
P (x ) =f (x 0) +f [x k , x k -1](x -x k ) +f [x k , x K -1, x k -2](x -x k )(x -x k -1)
有两个零点
x k +1=x k -
式中
2f (x k )
ω±
-4f (x k ) f [x k , x k _1, x k -2]
2
ω=f [x k , x k -1]+f [x k , x k -1, x k -2](x k -x k -1)
为了从式中定一个值x k +1，我们需要讨论根式前的正负号的取舍问题。
在x k , x k -1, x k -23个近似根中，自然假定x k 更接近所求的根x *，这时为了保证精度，我们选定上式中较接近x k 的一个值作为新的近似根x k -1，为此只要取根式前的符号与ω的符号相同即可。 抛物法的计算步骤
给定非线性方程f (x ) =0，误差界ε，迭代次数上限N ，则抛物法的计算步骤如下：
1） 计算ω=f [x k , x k -1]+f [x k , x k -1, x k -2](x k -x k -1) 2） 计算
2f (x k )
ω±ω-4f (x k ) f [x k , x k _1, x k -2]
x k +1=x k -
2
2
，代人
2f (x k )
ω±
-4f (x k ) f [x k , x k _1, x k -2]
得出的值后，再计算f (x k +1) 。
3） 若x k +1-x k ≤ε，则迭代停止，取x *≈x k +1；否则，令
(x k -2, x k -1, x k , f (x k -2), f (x k -1), f (x k )) =(x k -1, x k , x k +1, f (x k -), f (x k ), f (x k +1)) 4) 如果迭代次数k >N ，则认为该迭代格式对于所选的初值不收敛，迭代停止，否则重返步骤2）
下面是抛物法的MATLAB 函数代码： %抛物法求解非线性方程 function root=Parabola(f,a,b,x,eps) %f是非线性函数 %a为有根区间的下限 %b为有根区间的上限 %eps为根的精度 %root为求出的函数零点 %x为初始迭代点的值 if(nargin==4) eps=1.0e-4;
end
f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0) root=b; end if(f1*f2>0)
disp('两端点函数值乘积大于0！'); return; else tol=1;
fa=subs(sym(f),findsym(sym(f)),a); fb=subs(sym(f),findsym(sym(f)),b); fx=subs(sym(f),findsym(sym(f)),x); d1=(fb-fa)/(b-a); d2=(fx-fb)/(x-b); d3=(d2-d1)/(x-a); B=d2+d3*(x-b);
root=x-2*fx/(B+sign(B)*sqrt(B^2-4*fx*d3)); t=zeros(3); t(1)=a; t(2)=b; t(3)=x; while(tol>eps) t(1)=t(2); t(2)=t(3); t(3)=root;
f1=subs(sym(f),findsym(sym(f)),t(1)); f2=subs(sym(f),findsym(sym(f)),t(2)); f3=subs(sym(f),findsym(sym(f)),t(3)); d1=(f2-f1)/(t(2)-t(1)); d2=(f3-f2)/(t(3)-t(2)); d3=(d2-d1)/(t(3)-t(1)); B=d2+d3*(t(3)-t(2));
root=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3)); tol=abs(root-t(3)); end end 解：
在Matlab 命令窗口输入： root=Parabola ('f',-1,1,0.5) 回车运行得 ans =
0.56714
总 结
综合分析上述可得出各个算法的优缺点：
二分法是电子计算机上一种常用的算法，它的具有简单和易操作的优点，缺点是收敛较慢且不能求重根。牛顿迭代法具有至少平方收敛的速度，所以在迭代过程中只要迭代几次就会得到很精确的解，这是牛顿迭代法比简单迭代法优越的地方特别是当迭代点充分靠近精确解时，但选定的初值要接近方程的解，否则有可能得不到收敛的结果。再者, 牛顿迭代法计算量比较大，每次迭代要计算一次导数值，当表达式复杂，或无明显表达式时求解困难。对重根收敛速度较慢（线性收敛）。牛顿法是现在最常用的迭代方法。弦截法的收敛阶虽然低于Newton 法，但是迭代一次只要计算一次，不需要计算导数值，所以效率高，实际问题中经常使用。弦截法比牛顿迭代法收敛速度稍慢，但它的计算量比牛顿迭代法小，
特别当都函数的导数的计算比较复杂时，弦截法更显示了它的优越性。
综上所述，以上求解非线性方程的几种方法各有优缺点，通过Matlab 程序的实现可帮助我们更好理解它们的思想，在解决实际时，我们要根据实际情况选取适当的方法求解。
参考文献
[1] 李庆阳, 王能超，易大义. 数值分析. 第五版, 清华大学出版社，2008.12 [2] 张丰德.MATLAB 数值计算方法. 机械工业出版社, 2009. [3] 徐士良, 数值方法与计算机实现北京：清华大学出版社.2006年.