數(shù)值分析教教案18

Newton法1.算法的根本思想圖4-5牛頓法原理示意圖將非線性方程線性化，以線性方程的解逐步逼近非線性方程的解，這就是Newton法的根本思想。把函數(shù)在某一初始值點附近展開成Taylor級數(shù)有取其線性局部，近似地代替函數(shù)可得方程的近似式：設(shè)，解該近似方程可得：可作為方程〔4-1〕的近似解。重復(fù)以上過程，得迭代公式〔4-8〕按式〔4-8〕求方程〔4-1〕的近似解稱為Newton法。Newton法也是一種不動點迭代，其迭代函數(shù)為如圖〔4-5〕所示，從幾何上看，為曲線過點的切線，為切線與軸交點，那么是曲線上點處的切線與軸的交點。如此繼續(xù)下去，為曲線上點處的切線與軸的交點。因此Newton法是以曲線的切線與軸的交點作為曲線與軸的交點的近似，故Newton法又稱為切線法。2.切線法的收斂性理論可以證明，在有根區(qū)間上，連續(xù)且不變號，那么只要選取的初始近似根滿足，切線法必定收斂。它的收斂速度可如下推出。方程可以等價地寫成，假設(shè)移項可得。設(shè)，兩邊求導(dǎo)得，代入，那么必得。另一方面，比擬迭代公式和可知，。把函數(shù)在點展成泰勒級數(shù)，只取到二階導(dǎo)數(shù)那么有：由于，所以有，移項并注意，得出為了將式中的換成，對兩邊求導(dǎo)，并代入，那么有：將它代入前式得出〔4-9〕是個常數(shù)，式〔4-9〕說明用牛頓迭代公式在某次算得的誤差，與上次誤差的平方成正比，可見牛頓迭代公式的收斂速度很快，但計算實踐說明，當初值不夠好時，Newton法可能發(fā)散。一般可由問題的實際背景來預(yù)測或由對分區(qū)間法求得較好的初始值。3.Newton迭代公式Matlab實現(xiàn)按照算法4-4編寫迭代法的Matlab程序(函數(shù)名:Newton.m).function[p1,err,k,y]=newton(f,df,p0,delta,max1)%f是非線性函數(shù)%df是f的微商%p0是初始值%delta是給定允許誤差%max1是迭代的最大次數(shù)%p1是牛頓法求得的方程的近似值%err是p0的誤差估計%k是迭代次數(shù)%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所示圖形，即可知函數(shù)f(x)與x軸有交點，也就是說有根，并且從圖中能夠大致估算到根的位置。先用一個名為f.m的文件定義函數(shù)：functiony=f(x)y=x^3-3*x+2;再用一個名為df.m的文件定義函數(shù)的微商：functiony=df(x)y=3*x^2-3;然后在MATLAB命令窗口輸入：>>newton('f','df',1.2,10^(-6),10)圖4-6f(x)與x軸交點顯示圖回車得到如下結(jié)果〔只給出了最后兩次迭代結(jié)果〕：……k=9p1=1.0004err=4.1596e-004k=10p1=1.0002err=2.0802e-004ans=1.0002弦截法1.弦截法根本原理圖4-7弦截法原理示意圖用牛頓法解非性方程要知道的值，當比擬復(fù)雜時難以實現(xiàn)，下面要介紹的弦截法用在兩個點上的值構(gòu)造一次插值函數(shù)來回避微商的計算。給定非線性方程，選定曲線上的兩個點，，過這兩個點作一條直線，那么直線方程。當時，直線與軸交點為，這時用作為曲線與軸交點的近似值，顯然：這里為的精確解。然后用，，構(gòu)造直線，重復(fù)上述步驟，可以求出。如此進行下去，得到迭代格式：〔4-10〕迭代格式〔4-10〕實際上就是用差商取代牛頓公式中的微商的結(jié)果，所以弦截法可以看成牛頓法的一種變形。弦截法與牛頓雖然都是非線性方程線性化的方法，但牛頓法在計算時只用到的值，但是弦截法要用到和的值，故這就要給定兩個初值。切線法具有恒收斂且收斂速度快的優(yōu)點，但需要求出函數(shù)的導(dǎo)數(shù)；弦截法不需要求導(dǎo)數(shù)，收斂速度很快，但是需要知道兩個近似的初始值才能作出弦，要求的初始值條件較多。2.弦截法的MATLAB實現(xiàn)function[p1,err,k,y]=secant(f,p0,p1,delta,max1)%f是給定的非線性函數(shù)%p0,p1為初始值%delta為給定誤差界%max1是迭代次數(shù)的上限%p1為所求得的方程的近似解%err為p1-p0絕對值%k為所需要的迭代次數(shù)%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所示圖形，即可知函數(shù)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這就說明經(jīng)過了3次迭代得到滿足精足要求的近似解，且。圖4-8f(x)與x軸交點顯示圖拋物線法1.拋物法的根本原理設(shè)方程＝０的三個近似根,,，我們以這三點為節(jié)點構(gòu)造二次插值多項式的一個零點作為新的近似根，這樣確定的迭代過程稱拋物線法。在幾何圖形上，這種方法的根本思想是用拋物線與軸的交點作為所求根的近似位置〔圖6-7〕。圖4-9拋物法圖形現(xiàn)在推導(dǎo)拋物線法的計算公式．插值多項式有兩個零點：(4-11)式中為了從(4-11)式定出一個值，我們需要討論根式前正負號的取舍問題．在,,三個近似根中，自然假定更接近所求的根，這時，為了保證精度，我們選定(4-11)中較接近的一個值作為新的近似根，為此，只要取根式前的符號與ω的符號相同。2.拋物法的計算步驟給定非線性方程，誤差界ε，迭代次數(shù)上限N。〔1〕計算〔2〕計算，代入：得出的值后，計算?！?〕假設(shè)≤ε那么迭代停止，取≈；否那么的話，令：，n=n+1?！?〕如果迭代次數(shù)k＞N，那么認為該迭代格式對于所選初值不收斂，迭代停止；否那么的話重返步驟〔2〕。3.拋物法Matlab實現(xiàn)按照算法編寫迭代法的Matlab程序functionroot=Parabola(f,a,b,x,eps)%f是非線性函數(shù)%a為有根區(qū)間的左限%b為有根區(qū)間的右限%eps為根的精度%root為求出的函數(shù)零點%x為初始迭代點的值if(nargin==4)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);%subs命令是符號函數(shù)賦值f2=subs(sym(f),findsym(sym(f)),b);%findsym確定自由變量是對整個矩陣進行的if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于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】用拋物法求解方程在區(qū)間[1，4]內(nèi)的一個根。首先，在MA