数值分析教教案18

Newton法1.算法的根本思想图4-5牛顿法原理示意图将非线性方程线性化，以线性方程的解逐步逼近非线性方程的解，这就是Newton法的根本思想。把函数在某一初始值点附近展开成Taylor级数有取其线性局部，近似地代替函数可得方程的近似式：设，解该近似方程可得：可作为方程〔4-1〕的近似解。重复以上过程，得迭代公式〔4-8〕按式〔4-8〕求方程〔4-1〕的近似解称为Newton法。Newton法也是一种不动点迭代，其迭代函数为如图〔4-5〕所示，从几何上看，为曲线过点的切线，为切线与轴交点，那么是曲线上点处的切线与轴的交点。如此继续下去，为曲线上点处的切线与轴的交点。因此Newton法是以曲线的切线与轴的交点作为曲线与轴的交点的近似，故Newton法又称为切线法。2.切线法的收敛性理论可以证明，在有根区间上，连续且不变号，那么只要选取的初始近似根满足，切线法必定收敛。它的收敛速度可如下推出。方程可以等价地写成，假设移项可得。设，两边求导得，代入，那么必得。另一方面，比拟迭代公式和可知，。把函数在点展成泰勒级数，只取到二阶导数那么有：由于，所以有，移项并注意，得出为了将式中的换成，对两边求导，并代入，那么有：将它代入前式得出〔4-9〕是个常数，式〔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)fork=1:max1p1=p0-feval('f',p0)/feval('df',p0);err=abs(p1-p0);p0=p1;p1,err,k,y=feval('f',p1)if(err<delta)|(y==0),break,endend【例4-3】求方程的近似值，给定一个初始值，误差界为。首先，在MATLAB命令窗口输入：fplot('[x^3-3*x+2,0]',[-2.52.5]);grid;回车得到如图4-6所示图形，即可知函数f(x)与x轴有交点，也就是说有根，并且从图中能够大致估算到根的位置。先用一个名为f.m的文件定义函数：functiony=f(x)y=x^3-3*x+2;再用一个名为df.m的文件定义函数的微商：functiony=df(x)y=3*x^2-3;然后在MATLAB命令窗口输入：>>newton('f','df',1.2,10^(-6),10)图4-6f(x)与x轴交点显示图回车得到如下结果〔只给出了最后两次迭代结果〕：……k=9p1=1.0004err=4.1596e-004k=10p1=1.0002err=2.0802e-004ans=1.0002弦截法1.弦截法根本原理图4-7弦截法原理示意图用牛顿法解非性方程要知道的值，当比拟复杂时难以实现，下面要介绍的弦截法用在两个点上的值构造一次插值函数来回避微商的计算。给定非线性方程，选定曲线上的两个点，，过这两个点作一条直线，那么直线方程。当时，直线与轴交点为，这时用作为曲线与轴交点的近似值，显然：这里为的精确解。然后用，，构造直线，重复上述步骤，可以求出。如此进行下去，得到迭代格式：〔4-10〕迭代格式〔4-10〕实际上就是用差商取代牛顿公式中的微商的结果，所以弦截法可以看成牛顿法的一种变形。弦截法与牛顿虽然都是非线性方程线性化的方法，但牛顿法在计算时只用到的值，但是弦截法要用到和的值，故这就要给定两个初值。切线法具有恒收敛且收敛速度快的优点，但需要求出函数的导数；弦截法不需要求导数，收敛速度很快，但是需要知道两个近似的初始值才能作出弦，要求的初始值条件较多。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)fork=1:max1p2=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<delta)|(y==0),break,endend【例4-4】解非线性方程，给定初值为误差界为。首先，在MATLAB命令窗口输入：fplot('[x^3-x+2,0]',[-2.52.5]);grid;回车得到如图4-8所示图形，即可知函数f(x)与x轴有交点，也就是说有根，并且从图中能够大致估算到根的位置。解：先用一个名为f.m的文件定义；functiony=f(x)y=x^3-x+2;然后在MATLAB命令窗口输入：>>secant('f',-1.5,-1.52,10^(-6),11)k=0p0=-1.5000p1=-1.5200ans=0.1250ans=0.0082k=1p1=-1.5214err=0.0014y=-1.3633e-004k=2p1=-1.5214err=2.2961e-005y=1.4454e-007k=3p1=-1.5214err=2.4318e-008y=2.5460e-012ans=-1.5214这就说明经过了3次迭代得到满足精足要求的近似解，且。图4-8f(x)与x轴交点显示图抛物线法1.抛物法的根本原理设方程＝０的三个近似根,,，我们以这三点为节点构造二次插值多项式的一个零点作为新的近似根，这样确定的迭代过程称抛物线法。在几何图形上，这种方法的根本思想是用抛物线与轴的交点作为所求根的近似位置〔图6-7〕。图4-9抛物法图形现在推导抛物线法的计算公式．插值多项式有两个零点：(4-11)式中为了从(4-11)式定出一个值，我们需要讨论根式前正负号的取舍问题．在,,三个近似根中，自然假定更接近所求的根，这时，为了保证精度，我们选定(4-11)中较接近的一个值作为新的近似根，为此，只要取根式前的符号与ω的符号相同。2.抛物法的计算步骤给定非线性方程，误差界ε，迭代次数上限N。〔1〕计算〔2〕计算，代入：得出的值后，计算。〔3〕假设≤ε那么迭代停止，取≈；否那么的话，令：，n=n+1。〔4〕如果迭代次数k＞N，那么认为该迭代格式对于所选初值不收敛，迭代停止；否那么的话重返步骤〔2〕。3.抛物法Matlab实现按照算法编写迭代法的Matlab程序functionroot=Parabola(f,a,b,x,eps)%f是非线性函数%a为有根区间的左限%b为有根区间的右限%eps为根的精度%root为求出的函数零点%x为初始迭代点的值if(nargin==4)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);%subs命令是符号函数赋值f2=subs(sym(f),findsym(sym(f)),b);%findsym确定自由变量是对整个矩阵进行的if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('两端点函数值乘积大于0!');return;elsetol=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));endend【例4-5】用抛物法求解方程在区间[1，4]内的一个根。首先，在MA