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

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