第七講MATLAB中求方程的近似根(解)

1第七講MATLAB學習matlab準解析法、數(shù)值方法以及迭代方法，把握對分法、迭代法、牛頓切法線求方程近似根的根本過程；把握求代數(shù)方程〔組〕的解的求解命令．教學重點：求方程近似解的幾種迭代方法，代數(shù)方程〔組〕的解的求解命令的使用程〔組〕的求解命令及使用技巧．教學難點：方程的近似求解和非線性方程〔組〕的求解．一、問題背景和試驗?zāi)康膄(x)0的根是最常見的數(shù)學問題之一〔這里稱為代數(shù)方程，主要是想和后面的微分方程區(qū)分開．為簡明起見，在本試驗的以下表達中，把代數(shù)方程簡稱為方程f(x)f(x)0為線性方程，否則稱之為非線性方程．f(x)0f(x)的多樣性，尚無一般的解析解法可使用，但假設(shè)能滿足實際要求．同時對于多未知量非線性方程(組)而言，簡潔的迭代法也是可以做出來的，但在這里我們介紹相關(guān)的命令來求解，不用迭代方法求解．通過本試驗，到達下面目的：了解對分法、迭代法、牛頓切線法求方程近似根的根本過程；求代數(shù)方程〔組〕的解．首先，我們先介紹幾種近似求根有關(guān)的方法：對分法到滿足精度為止．對分法適用于求有根區(qū)間內(nèi)的單實根或奇重實根．f(x)在[a,b上連續(xù)，f(af(b)0，即f(a)0，f(b)0f(a)0，f(b)0依據(jù)連續(xù)函數(shù)的介值定理，在(ab內(nèi)至少存在一點f()0．下面的方法可以求出該根：x0

(ab/2f(x；0f(x0

)0x0

f(x)0xx．0假設(shè)f(af(x0

)0a1

a，b1

xf(af(x0

)0a1

x，b0

b；x(a1

b/2．……，有a1

、bxk

(ak

b)/2．kf(xk

) 為預(yù)先給定的精度要求)，退出計算，輸出結(jié)果xk

(ak

b)/2；k反之，返回(1)，重復(fù)(1)，(2)，(3)．以上方法可得到每次縮小一半的區(qū)間序列{[abk k

]}，在(abk k

中含有方程的根．當區(qū)間長bk

axk

(ak

b/2為根的近似值，明顯有kx (bk k

a)/2(bk

a

)/(22)

(ba)/2k1以上公式可用于估量對分次數(shù)k．分析以上過程不難知道，對分法的收斂速度與公比為

1的等比級數(shù)一樣．由于22101024，可知大約對分10它常用來摸索實根的分布區(qū)間，或求根的近似值．迭代法f(x)0構(gòu)造一個等價方程x(x)．則迭代方程是：x

(1)x(x)，1/(1”(x

))，其中”(x)1．

k k k k k k斂．Altken松弛法要先計算”(xk

，在使用中有時不便利，為此進展出以下的Altken公式：x(1)k

(x) ；kx(2)k

(x(1))；kx x(2)(x(2)x(1))2/(x(2)

2x(1)x

) k1 k k

k k kAltken5)．牛頓(Newton)法〔牛頓切線法〕f(x)0是非線性方程其迭代公式為：

x x(f(x)/f”(x)) k0,1,2,k k k牛頓法的缺點是：(1)對重根收斂很慢；(2)對初值x要求較嚴，要求x相當接近真值0 0x*．因此，常用其他方法確定初值x，再用牛頓法提高精度．0以下是本試驗中的幾個具體的試驗1：對分法x33x10的實根的分布區(qū)間，再利用對分法在這些區(qū)間上分別求出根的近似值．程序如下：function[y,p]=erfenclc,x=[];a=[];b=[];a(1)=1;b(1)=2;i=1;x(i)=(a(i)+b(i))/2;e=abs(f(x(i)));ezplot(”x^3-3*x+1”,[a(1),b(1)]);holdon, whilee>10^(-5)plot([a(i),a(i)],[0,100],[x(i)x(i)],[0100],[b(i)b(i)],[0100]),pause(0.5)iff(a(i))*f(x(i))<0a(i+1)=a(i);b(i+1)=x(i);x(i+1)=(a(i+1)+b(i+1))/2;elseend

a(i+1)=x(i);b(i+1)=b(i);x(i+1)=(a(i+1)+b(i+1))/2;e=abs(f(x(i)));i=i+1;endy=x(i);p=[a;x;b]”functionu=f(x)u=x^3-3*x+1;endend圖形如下:1.532132：一般迭代法承受迭代過程：x (x)求方程x33x10在0.5四周的根，準確到第4位小k數(shù)．x(x31)/3用迭代公式：x (x31)/3，k0,1,2,k具體試驗3：迭代法的加速1——松弛迭代法迭代公式為

(x)(x31)/3，”(x)x2，k

1/(1x2)kxk1

(1)xk k

(xk

31)/3clc;x=[];w=[];x(1)=1;w(1)=1/(1-x(1));fori=1:10w(i)=1/(1-x(i)); x(i+1)=(1-w(i))*x(i)+w(i)*(x(i)^3+1)/3;endx另外有程序可以參考，詳見參見附錄4．具體試驗4：迭代法的加速2——Altken迭代法迭代公式為：

x(1)

(x31)/3，x(2)

(x(1)31)/3k k k kx x(2)(x(2)x(1))2/(x(2)

2x(1)x

〕symsxfxgx;

k1 k k k

k k kgx=(x^3+1)/3;fx=x^3-3*x+1;disp(”k x x1 x=0.5;k=0;ffx=subs(fx,”x”,x);whileabs(ffx)>0.0001;u=subs(gx,”x”,x);v=subs(gx,”x”,u);disp([num2str(k),” ”,num2str(x),” ”,num2str(u),” x=v-(v-u)^2/(v-2*u+x);k=k+1;ffx=subs(fx,”x”,x);enddisp([num2str(k),” ”,num2str(x),” ”,num2str(u),” ”,num2str(v)])%〔數(shù)值計算〕function[y,p]=althken%求方根的迭代程序x(1)=6;i=1;p=[];ezplot(”x^3-3*x+1”,[x(1)-9,x(1)+1]);holdonplot([x(1)-20,x(1)+2],[0,0])whileabs(f(x(i)))>=10^(-5)plot(x(i),0,”*”)4end

t1=phi(x(i));t2=phi(t1); x(i+1)=t2-(t2-t1)^2/(t2-2*t1+x(i)+eps);p=[p;[i,x(i),t1,t2]]; i=i+1; pause(0.1)p,y=x(i),i,formatfunctionu=phi(x)u=(x^3+1)/3;endfunctionu=f(x)u=x^3+1-3*x;endend5：牛頓法x33x10在-22f(x)x33x1f”(x)3x23迭代公式：function[y,p]=newton%求方根的迭代程序

x

x2(xk

33xk

1)/(3x23)ki=1;p=[];ezplot(”x^3-3*x+1”,[x(1)-9,x(1)+1]);holdonplot([x(1)-20,x(1)+2],[0,0])whileabs(f(x(i)))>=10^(-5)plot(x(i),0,”*”), p=[p;[i,x(i)]]; i=i+1; pause(0.1)endformatshort,p,y=x(i),i,functionu=df(x)u=3*x^2-3;endfunctionu=f(x)u=x^3+1-3*x;endend結(jié)果：5結(jié)果為：1.5321※進一步思考：用迭代法求3的平方根.迭代公式為x (x3/x)/2.編寫M函數(shù)n1 n nMy_sqrt.m,3x.要求誤差小于105．僅要求寫出源程序．試使用以上介紹的迭代法來相互比較參考程序：functiony=my_sqrt(a)%求方根的迭代程序ifnargin~=1|~isa(a,”double”), 輸入數(shù)字為一個正數(shù)!”),endifa<0, error(”輸入數(shù)字為正數(shù)!”),endifa>0formatlonge, x(1)=0; x(2)=1;i=1;whileabs(x(i+1)-x(i))>=10^(-5)i=i+1;x(i+1)=1/2*(x(i)+a/(x(i)+eps));endy=x(i+1);i,formatend

現(xiàn)在我們簡潔介紹圖解法如何來求解一元方程和二元方程的根：例：exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)=0.5>>ezplot(”exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5”,[05])>>holdon,line([0,5],[0,0])t=3.5203>>symsx;t=3.5203;vpa(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans=610-.43167073997540938989914138801396e-4x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)=0y^2*cos(y+x^2)+x^2*exp(x+y)=0>>ezplot(”x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)”)>>holdonezplot(”y^2*cos(y+x^2)+x^2*exp(x+y)”)具體的結(jié)果請大家自己下來運行〔命令〕求解方程及簡介solve(”f(x)”)，f(x)為一個具體的表達式．roots(A)，Axfzero(”f(x)”,x0)，f(x)為一個具體的表達式，x0linsolve(A,b)，A，b單變量的多項式方程求根：命令格式：roots(A)例：x^3-6*(x^2)-72*x-27=0;>>p=[1-6-72-27]>>r=roots(p)r=12.1229-5.7345-0.3884多項式型方程的準解析解法命令格式：[x,…]=solve(eqn1,eqn2,…)例：x^2+y^2-1=00.75*x^3-y+0.9=0>>symsxy;>>[x,y]=solve(”x^2+y^2-1=0”,”75*x^3/100-y+9/10=0”)檢驗：>>[eval(”x.^2+y.^2-1”),eval(”75*x.^3/100-y+9/10”)]具體結(jié)果就請大家下來自己運行線性方程組的求解例：求線性方程組mxb的解，m=[12345;23456;345678;45678;56780]，b=[1;2;3;4;5]fori=1:5forj=1:5m(i,j)=i+j-1;endendm(5,5)=0;b=[1:5]”;linsolve(m,b)非線性方程數(shù)值求解用格式為：z=fzero(”fname”,x0,tol,trace)其中fnamex0toltol=eps，trace指10trace=0．例:f(x)=x-10x+2=0x0=0.5步驟如下：functionfx=funx(x)fx=x-10.^x+2;z=fzero(”funx”,0.5)z= 0.3758非線性方程組的求解．fsolveX=fsolve(”fun”,X0,option)X020戶可以使用optimset命令將它們顯示出來．假設(shè)想轉(zhuǎn)變其中某個選項，則可以調(diào)用optimset函數(shù)來完成．例如，Display‘off’‘iter’‘final’optimset(‘Display’,‘off’)Display‘off’．例:求以下非線性方程組在(0.5,0.5)functionq=myfun(p)x=p(1);y=p(2);q(1)=x-0.6*sin(x)-0.3*cos(y);q(2)=y-0.6*cos(x)+0.3*sin(y);fsolvex=fsolve(”myfun”,[0.5,0.5]”,optimset(”Display”,”off”))x= 0.63540.3734q=myfun(x)q= 1.0e-009*0.2375 0.2957可見得到了較高精度的結(jié)果．精品案例：螺旋線與平面的交點問題 ：螺旋線與平面相交的狀況多種多樣,依據(jù)螺旋線與平面方程的不同可以相交,也可以不相交.在相交的狀況下,可以交于一點,也可以交于好多點.對于各種相交的狀況,要求其交點的坐標并不是一件簡潔的事.本次試驗就以此為背景爭論下面的具體問題：螺旋線的參數(shù)方程為x4cos,y4sin,z,08.xy0.5z20求該螺旋線與平面的交點.要求：求出全部交點的坐標；在同一圖形窗口畫出螺旋線、平面和交點.試驗過程：問題分析fsolve等.先對方程化簡,削減變量個數(shù),使用圖解方法求方程的根.再分別畫出螺旋線,平面,及其交點.算法描述與分析fsolve,選擇適當?shù)某踔?求其數(shù)值解;再分別會出圖形；最終對圖形作出必要的修飾.源程序及注釋cos4sin0.520cos4sin20.5建立下面M文件intersect_point.m%使用圖解法求交點，并且三維圖%畫圖確定解的個數(shù)和或許位置theta=0:0.01:8*pi;y1=4*(cos(theta)+sin(theta));y2=2-0.5*theta;plot(theta,y1,theta,y2)%畫出兩個函數(shù)的圖形%畫螺旋線%theta=0:pi/100:8*pi;x=4*cos(theta);y=4*sin(theta);z=theta;figure %建圖形窗口plot3(x,y,z) holdon%透亮的畫平面%x1=-5:0.1:5;[-4,4]有關(guān).y1=x1;[X1Y1]=meshgrid(x1,y1);%網(wǎng)格化，畫曲面Z1=4-2*X1-2*Y1;surf(X1,Y1,Z1) mesh(X1,Y1,Z1)shadingflatalpha(0.5) %設(shè)置透亮度alpha(”z”) %設(shè)置透亮度方向ttheta%i=1forn=[2,5,9,11]%依據(jù)畫圖確定解的或許位置作為初值t(i)=fsolve(inline(”4*cos(t)+4*sin(t)+0.5*t-2”),n)%選擇不同初值求交點x0(i)=4*cos(t(i));y0(i)=4*sin(t(i));z0(i)=t(i);i=i+1;endplot3(x0,y0,z0,”ro”)測試結(jié)果〔寫清輸入輸出狀況〕交點坐標為：

08

內(nèi)與三角曲線有4個交點.theta的數(shù)值解為：t=[2.1961 5.3759 9.1078 11.1023]四個交點的近似坐標為：x0=[-2.3413 2.4635 -3.8007 0.4261]y0=[3.2432-3.15141.2468-3.9772]z0=[2.19615.37599.107811.1023]調(diào)試和運行程序過程中產(chǎn)生的問題及實行的措施求交點的時候會消滅重根和漏根的情形,通過選擇適當?shù)某踔当荛_了上述狀況.對算法和程序的爭論、分析,改進設(shè)想及其它閱歷教訓(xùn)solve函數(shù)只能求解一個數(shù)值解,不能全部求出;fsolve函數(shù)好;為了滿足更好的視覺1112效果,可以對圖形進展進一步的修飾.習題1.f(x)3x5x42x3x23sin(xyyex0(1)x2y2 (2)求解方程:cos(x)xex求解多項式方程x9810求以下代數(shù)方程〔組〕的解：(1) x5x10(2) 2x3y0 ①4x23y1 ②〔1一般迭代法〔2與之相應(yīng)的松弛迭代法和x33x10在1.4四周的根，準確到4變化．Altken迭代法、牛頓切法線等5方程txsin(x)的正的近似根，0t1．(建議取t0.5．時間許可的話，可進一步考慮t0.25的狀況．)五、附錄為供近似求根的算法1：對分法程序〔fulu1.m〕symsxfx;a=0;b=1;fx=x^3-3*x+1;x=(a+b)/2;k=0;ffx=subs(fx,”x”,x);ifffx==0;disp([”therootis:”,num2str(x)])elsedisp(”kakbkf(xk)”)whileabs(ffx)>0.0001&a<b;num2str(a),””,num2str(b),””,num2str(ffx)])fa=subs(fx,”x”,a);ffx=subs(fx,”x”,x);iffa*ffx<0b=x;elsea=x;endk=k+1;x=(a+b)/2;enddisp([num2str(k),” ”,num2str(a),” ”,num2str(b),” ”,num2str(ffx)])end注：試驗時，可將第2行的a、b改為其它區(qū)間端點進展其它試驗．2：一般迭代法〔fulu2.m〕symsxfxgx;gx=(x^3+1)/3;fx=x^3-3*x+1;disp(”k x x=0.5;k=0;ffx=subs(fx,”x”,x);whileabs(ffx)>0.0001;disp([num2str(k),” ”,num2str(x),” ”,num2str(ffx)]);x=subs(gx,”x”,x);ffx=subs(fx,”x”,x);k=k+1;enddisp([num2str(k),” ”,num2str(x),” ”,num2str(ffx)])3：收斂/發(fā)散推斷〔fulu3.m〕symsxg1g2g3dg1dg2dg3;x1=0.347;x2=1.53;x3=-1.88;g1=(x^3+1)/3;dg1=diff(g1,”x”);g2=1/(3-x^2);dg2=diff(g2,”x”);g3=(3*x-1)^(1/3);dg3=diff(g3,”x”);disp([”1 ”,num2str(abs(subs(dg1,”x”,x1))),” ”,...num2str(abs(subs(dg1,”x”,x2))),” ”,num2str(abs(subs(dg1,”x”,x3)))])disp([”2 ”,num2str(abs(subs(dg2,”x”,x1))),” ”,...num2str(abs(subs(dg2,”x”,x2))),” ”,num2str(abs(subs(dg2,”x”,x3)))])disp([”3 ”,num2str(abs(subs(dg3,”x”,x1))),” ”,...num2str(abs(subs(dg3,”x”,x2))),” ”,num2str(abs(subs(dg3,”x”,x3)))])4：松弛迭代法〔fulu4.m〕symsfxgxxdgx;gx=(x^3+1)/3;fx=x^3-3*x+1;dgx=diff(gx,”x”);x=0.5;k=0;ggx=subs(gx,”x”,x);ffx=subs(fx,”x”,x);dgxx=subs(dgx,”x”,x);disp(”k x w”)whileabs(ffx)>0.0001;w=1/(1-dgxx); disp([num2str(k),” ”,num2str(x),” ”,num2str(w)])x=(1-w)*x+w*ggx;k=k+1;ggx=subs(gx,”x”,x);ffx=subs(fx,”x”,x);dgxx=subs(dgx,”x”,x);enddisp([num2str(k),” ”,num2str(x),” ”,nu