第七講MATLAB中求方程的近似根(解)_第1頁
第七講MATLAB中求方程的近似根(解)_第2頁
第七講MATLAB中求方程的近似根(解)_第3頁
第七講MATLAB中求方程的近似根(解)_第4頁
第七講MATLAB中求方程的近似根(解)_第5頁
已閱讀5頁,還剩12頁未讀, 繼續(xù)免費閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進行舉報或認領(lǐng)

文檔簡介

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è)想及其它閱歷教訓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

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負責。
  • 6. 下載文件中如有侵權(quán)或不適當內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論