数值分析教教案18_数值分析教教案

数值分析教教案18由刀豆文库小编整理，希望给你工作、学习、生活带来方便，猜你可能喜欢“数值分析教教案”。

4.1.3 Newton法

1.算法的基本思想

图4-5 牛顿法原理示意图

将非线性方程线性化，以线性方程的解逐步逼近非线性方程的解，这就是Newton法的基本思想。

把函数f(x)在某一初始值x0点附近展开成Taylor级数有

f(x0)f(x)f(x0)(xx0)f(x0)(xx0)取其

2!线性部分，近似地代替函数f(x)可得方程的近似式：

2f(x)f(x0)(xx0)f(x0)0

设f(x0)0，解该近似方程可得：

f(x0)x1x0f(x0)

可作为方程（4-1）的近似解。重复以上过程，得迭代公式

xk1f(xk)xk(k0,1,2,)（4-8）

f(xk)按式（4-8）求方程（4-1）的近似解称为Newton法。Newton法也是一种不动点迭代，其迭代函数为

f(x)g(x)xf(x)

如图（4-5）所示，从几何上看，yf(x0)f(x0)(xx0)为曲线yf(x)过点(x0,f(x0))的切线，x1为切线与x轴交点，x2则

x轴的交点。如此继续下去，xk1为曲线上点(x,f(x))处的切线与x轴的交点。因此Newton法是以曲线的切线与x轴的交点作为曲线与x轴的交点的近似，故是曲线上点(x1,f(x1))处的切线与

kkNewton法又称为切线法。

2.切线法的收敛性

理论可以证明，在有根区间[a,b]上，f(x)0,f(x)连续且不变号，则只要选取的初始近似根

x0满足f(x0)f(x0)0，切线法必定收敛。它的收敛速度可如下推出。

方程f(x)0可以等价地写成f(x)(xx)f(x)，若f(x)0f(x)f(x)移项可得xxf(x)。设g(x)xf(x)，两边求导得f(x)f(x)g(x)f(x)0g(x)0。xx，则必得 2，代入[f(x)]另一方面，比较迭代公式xk1f(xk)xkf(xk)和f(x)g(x)xg(x)xg(x)x可知。把函数在点展k1kf(x)成泰勒级数，只取到二阶导数则有：

g(x)2g(xk)g(x)g(x)(xkx)(xkx)由

2g(x)g(xk)g(x)(xkx)2，移

2xk1g(x)0，所以有xk1于xg(x)，得出 项并注意

g(x)xk1g(x)xk1x(xkx)2

2为了将式中的g(x)换成ff(x)f(x)(x)，对g(x)[f(x)]2两边求导，并代入f(x)0，则有：

f(x)g(x)f(x)*将它代入前式得出

f(x)2xk1x(xkx)（4-9）2f(x)f(x)2f(x)是个常数，式（4-9）表明用牛顿迭代公式在某次算得的误差，与上次误差的平方成正比，可见牛顿迭代公式的收敛速度很快，但计算实践表明，当初值不够好时，Newton法可能发散。一般可由问题的实际背景来预测或由对分区间法求得较好的初始值。

3.Newton迭代公式Matlab实现

按照算法4-4编写迭代法的Matlab程序(函数名:Newton.m).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

3x3x2的近似值，给定一个初始值【例4-3】求方程p01.2，误差界为106。

首先，在MATLAB命令窗口输入：

fplot('[x^3-3*x+2,0]',[-2.5 2.5]);grid;回车得到如图4-6所示图形，即可知函数f(x)与x轴有交点，也就是说有根，并且从图中能够大致估算到根的位置。

3f(x)x3x2： 先用一个名为f.m的文件定义函数function y=f(x)y=x^3-3*x+2;

2df(x)3x3： 再用一个名为df.m的文件定义函数的微商function y=df(x)y=3*x^2-3;然后在MATLAB命令窗口输入： >> newton('f','df',1.2,10^(-6),10)

图4-6 f(x)与x轴交点显示图

回车得到如下结果（只给出了最后两次迭代结果）： …… k = 9 p1 = 1.0004 err = 4.1596e-004 k = 10 p1 = 1.0002 err = 2.0802e-004 ans = 1.0002 4.1.4弦截法

1.弦截法基本原理

图4-7弦截法原理示意图

f(xk)的值，当f比较复杂时难以实现，下面要介绍的弦截法用f(x)在两个点上的值构造一次插值函用牛顿法解非性方程要知道数来回避微商的计算。

给定非线性方程f(x)0，选定曲线

yf(x)上的两个点P0(x0,f(x0))，P1(x1,f(x1))，过这两个点作一条直线P0P1，则直f(x1)f(x0)(xx1)。当线方程yf(x1)x1x0f(x0)f(x1)时，f(x1)(x1x0)x2x1xx直线P0P与轴交点为，这时用2作为1f(x1)f(x0)曲线yf(x)与x轴交点的近似值，显然：

x2xmin(x1x,x0x)

这里***x*为f(x)0的精确解。然后用

P1(x1,f(x1))，P2(x2,f(x2))，构造直线P1P2，重复上述步骤，可以求出x3。

如此进行下去，得到迭代格式：

f(xk)xk1xk(xkxk1)(k0,1,)（4-10）f(xk)f(xk1)f(xk)f(xk1)迭代格式（4-10）实际上就是用差商取代牛顿公xkxk1式xk1f(xk)xk(x)的结果，所以弦截法可以看成f中的微商f(xk)牛顿法的一种变形。

弦截法与牛顿虽然都是非线性方程线性化的方法，但牛顿法在计算xk1时只用到xk的值，但是弦截法要用到xk和xk1的值，故这就要给定两个初值。切线法具有恒收敛且收敛速度快的优点，但需要求出函数的导数；弦截法不需要求导数，收敛速度很快，但是需要知道两个近似的初始值才能作出弦，要求的初始值条件较多。2.弦截法的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

xx20，给定初值为

3p01.5,p11.52,误差界为106。

首先，在MATLAB命令窗口输入：

fplot('[x^3-x+2,0]',[-2.5 2.5]);grid;回车得到如图4-8所示图形，即可知函数f(x)与x轴有交点，也就是说有根，并且从图中能够大致估算到根的位置。

3f(x)xx2； 解：先用一个名为f.m的文件定义function y=f(x)y=x^3-x+2;然后在MATLAB命令窗口输入：

>> secant('f',-1.5,-1.52,10^(-6),11)k =0 p0 =-1.5000 p1 =-1.5200 ans =0.1250 ans =0.0082 k =1 p1 =-1.5214 err =0.0014 y =-1.3633e-004 k =2 p1 =-1.5214 err =2.2961e-005 y =1.4454e-007 k =3 p1 =-1.5214 err =2.4318e-008 y =2.5460e-012 ans =-1.5214 这就表明经过了3次迭代得到满足精足要求的近似解xx31.5214，且f(x3)2.5461012。

图4-8 f(x)与x轴交点显示图

4.1.5抛物线法

1.抛物法的基本原理

设已知方程f(x)＝０的三个近似根xk,xk-1,xk-2，我们以这三点为节点构造二次插值多项式P(x)的一个零点xk+1作为新的近似根，这样确定的迭代过程称抛物线法。

在几何图形上，这种方法的基本思想是用抛物线与

x轴的交点xk+1作为所求根的近似位置（图6-7）。

图4-9抛物法图形

现在推导抛物线法的计算公式．插值多项式

P(x)=f(xk)+ f[xk,x k-1](x-xk)+f[xk,xk-1,xk-2](x-xk)(x-xk-1)

有两个零点：

xk+1=xk－式中

2f(xk)ω±ω2-4f(xk)f[xk,xk-1,xk-2](4-11)

ω=f[xk,xk-1]+f[xk,xk-1,xk-2](xk-xk-1)

为了从(4-11)式定出一个值取舍问题． 在xk+1，我们需要讨论根式前正负号的xk,xk－1,xk－2三个近似根中，自然假定xk更接近所求的根x*，这时，为了保证精度，我们选定(4-11)中较接近xk的一个值作为新的近似根xk+1，为此，只要取根式前的符号与ω的符号相同。

2.抛物法的计算步骤 给定非线性方程f(x)0，误差界ε，迭代次数上限N。（1）计算ω=（2）计算f[xk,xk－ 1]+f[xk,xk－ 1,xk－ 2](xk－xk－ 1)

2f(xk)，代入：

ω±ω2-4f(xk)f[xk,xk-1,xk-2]xk+ 1=xk-2f(xk)ω±ω2-4f(xk)f[xk,xk-1,xk-2]

得出xk+1的值后，计算f(xk+1)。

*x－xx（3）若k+1≈xk+1；否则的话，令： k≤ε则迭代停止，取(xk－ 2,xk－ 1,xk,f(xk－ 2),f(xk－ 1),f(xk))=(xk－ 1,xk,xk+ 1,f(xk－ 1),f(xk),f(xk+ 1))，n=n+1。

（4）如果迭代次数k＞N，则认为该迭代格式对于所选初值不收敛，迭代停止；否则的话重返步骤（2）。3.抛物法Matlab实现

按照算法编写迭代法的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);% subs命令是符号函数赋值

f2=subs(sym(f),findsym(sym(f)),b);% findsym确定自由变量是对整个矩阵进行的 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 【例4-5】用抛物法求解方程lgx+x=2在区间[1，4]内的一个根。

首先，在MATLAB命令窗口输入：

fplot('[sqrt(x)+log(x)-2]',[1 4]);grid;回车得到如图4-10所示图形，即可知函数f(x)与x轴有交点，也就是说有根，并且从图中能够大致估算到根的位置。

1.510.50-0.5-1

图4-10 f(x)与x轴交点显示图

然后在MATLAB命令窗口输入：

>> r=Parabola(' sqrt(x)+log(x)-2',1,4,2)回车得到如下结果： r = 1.8773 从结果中可以知道，方程lgx+11.522.533.54x=2的一个根x=1.8773。