数值计算课程设计报告非线性方程求根文档格式.docx
《数值计算课程设计报告非线性方程求根文档格式.docx》由会员分享，可在线阅读，更多相关《数值计算课程设计报告非线性方程求根文档格式.docx（26页珍藏版）》请在冰豆网上搜索。
三、基本算法思想与实现
（1）二分法
单变量函数方程：
f（x）=0
其中，f（x）在闭区间[a，b]上连续、单调，且f（a）*f（b）<
0,则有函数的介值定理可知，方程f（x）=0在（a，b）区间有且只有一个解
，二分法是通过函数在区间端点的符号来确定
所在区域，将有根区间缩小到充分小，从而可以求出满足给定精度的根
的近似值。
下面研究二分法的几何意义：
设
=1,
=b,区间
，中点
=
及
，若
=0，则
，若f（
）*f（
）<
0，令
，
，则根
[
]中，这样就得到长度缩小一半的有根区间[
]，若f（
],即f（
）f（
0,此时
-
对有根区间[
]重复上述步骤，即分半求中点，判断中电处符号，则可得长度有缩小一半的有根区间[
],
如图所示：
重复上述过程，第n步就得到根
的近似序列
及包含
的区间套，如下：
（1）
（2）
（3）
=…=
（4）
且|
|
（n=1,2,3…..）
显然lim
，且
以等比数列的收敛速度收敛于
，因此用二分法求f（x）=0的实根
可以达到任意指定精度。
（2）牛顿迭代法
设方程f（x）=0在其根
的某个领域U（
）有一阶连续导数，且f’（
）≠0。
求f（x）=0的根
，首先要将f（x）=0转化为等价形式
，并使
（x）满足不动点迭代的一般理论。
于是我们令
（x）=x+h（x）f（x）,可由
‘（
）=0来确定h（x）的结构，根据
’（x）=1+h’（
）+h（
）f’（x1）=1+h（
）f’（
）=0可得
h（
）=-1/f’（
），由于f’（x）≠0，且f’（x）连续，因此当h（x）=-1/f’（x）时，h’（x1）=0，即令
（x）=x-f（x）/f‘（x），从而有迭代格式
=
（k=0,1,2，…..）
由于
…….都在U领域里，从而当B比较小时，可用f’（
）可近似代替f’（
）,
=
-
，此方法称为牛顿迭代法
下面研究牛顿法的几何意义：
设r是方程f（x）=0的根，选取
作为的r初始近似值，经过（
f（
））做曲线y=f（x）的切线的方程：
y=f（
）+f’（
）（x-
），求出L与x的交点的横坐标
-f（
）/f’（
）,称
为r的一次近似值经过点（
））做切线y=f（x）的切线，并求出该切线与x轴的交点横坐标：
称为r的二次近似值，重复以上操作可以得到r的近似值序列。
下述三个定理分别讨论了牛顿法的收敛性质：
定理1：
对于方程f（x）=0，设f（x）在[a，b]上有二阶连续导数且满足下述条件：
（1）f（a）f（b）<
0;
（2）f’（x）
0,
0,对任意的x
[a,b];
（3）存在
[a,b]，使f（
）
>
0，
则由牛顿法产生的迭代序列
收敛于f（x）=0的根
定理2：
（2）对任意的x
[a,b]，f’（x）
<
b-a,
b-a
则对于任何
[a,b]，由牛顿法产生的迭代序列
定理3：
是方程f（x）=0的根，在
的某个开区间
连续且f’（x）
0，则存在
0,当
【
-,
+
】时，由牛顿迭代法
（k=0,1,2，…..）式产生的序列
是以不低于二阶的收敛速度收敛到
.
（3）弦截法
为方程f（x）=0的两个近似根。
用差商得：
f（
）-f（
）/
代替牛顿迭代公式中的导数f’（
），于是得到如下的迭代公式：
。
下面研究割弦法的几何意义：
经过点（
））及点（
））两点作割线，其点斜式方程为：
Y=f（
）-
，其零点为X=
把X用
表示即得到迭代格式，它又称为双点弦割法，需要两个初值
此割线与X轴交点的横坐标就是新的近似值
，所以弦截法又称为割线法，如图所示。
下面三个定理为弦割法收敛定理：
设f（x）在其零点
的邻域U（
）=[
]（
0）有二阶连续导数，
，则当
U（
）时，由割弦法式产生的序列
收敛于
，且收敛的阶为1.618。
在区间[a,b]上连续，且满足下述三点
[a,b]，有f’（x）
则对于任意初始
[a,b]，由弦割法产生的迭代序列
收敛于f（x）=0唯一的根
设在其零点
）=[
]（
0）有二阶连续导数，f’（x）
0则当
U（
）时，由弦割
式产生的序列
四、具体应用实例分析
求解
在
附近的根。
建立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为函数
ifya*yb>
disp（'
注意：
ya*yb>
0,请重新调整区间端点a和b.'
）,
return
end
max1=-1+ceil（（log（b-a）-log（abtol））/log
（2））;
%ceil是向
方向取整
fork=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]
Ifyx==0
a=x;
b=x;
else
ifyb*yx>
yb=yx;
ya=yx;
ifb-a<
abtol,
return,
k=max1;
x;
wuca;
建立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
建立newtonqx-M文件：
function[k,xk,yk,piancha,xdpiancha]=newtonqx（x0,tol,ftol,gxmax）
x
（1）=x0;
fori=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;
[（i-1）xkykpianchaxdpiancha]
return;
ifi>
gxmax
请注意：
迭代次数超过给定的最大值gxmax。
'
[（i-1）,xk,yk,piancha,xdpiancha]'
;
建立FNQ原函数文件：
functiony=fnq（x）
建立DFNQ导函数文件：
functiony=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；
升级会员 