數(shù)學(xué)實(shí)驗(yàn)之2.方程求解(1)(2004.3)

1、 方程求解方程求解(1)(1)上一頁下一頁主 頁1 1 復(fù)習(xí)求解方程的基本原理和方法；復(fù)習(xí)求解方程的基本原理和方法；2 2 掌握迭代算法；掌握迭代算法；3 3 熟悉熟悉MATLABMATLAB軟件編程環(huán)境；掌握軟件編程環(huán)境；掌握MATLABMATLAB編程語句；編程語句；4 4 了解迭代過程的圖形表示，分形與混沌學(xué)科等；了解迭代過程的圖形表示，分形與混沌學(xué)科等；5 5 通過范例展現(xiàn)求解實(shí)際問題的初步建模過程；通過范例展現(xiàn)求解實(shí)際問題的初步建模過程；返 回主要內(nèi)容主要內(nèi)容方程的常用求解方法方程的常用求解方法迭代過程的圖解迭代過程的圖解實(shí)實(shí) 驗(yàn)內(nèi)容驗(yàn)內(nèi)容引言及引例引言及引例圖形放大簡(jiǎn)單迭代牛頓迭代

2、上一頁下一頁主 頁 求高次代數(shù)方程的根是一個(gè)基本的、古老的數(shù)學(xué)問題。0011102nnnnaxaxaxacbxax高階方程能否求解？定義：一些由實(shí)際問題列出的方程中常常包含三角函數(shù)、指數(shù)函數(shù)等，它與n(2)代數(shù)方程一起統(tǒng)稱為非線性方程非線性方程(組組)，記作0)(xf引言例：exsin(x)+5x3=0一元二次方程一元二次方程高階多項(xiàng)式方程高階多項(xiàng)式方程上一頁下一頁主 頁 計(jì)算機(jī)技術(shù)和數(shù)學(xué)軟件的飛速發(fā)展，為一些古老數(shù)學(xué)問題的解決提供了極大的方便。如利用MATLAB的圖形功能就能幫助我們判斷方程有沒有根，確定根的近似位置，例如：0e50305813261282345xxxxxxx-1-0.500

3、.511.522.53-60-40-20020406080100120140 x8 x5-12 x4-26 x3-13 x2+58 x+308x5-12x4-26x3-13x2+58x+30 00.511.522.53-4-3-2-10123x5 x-exp(x)5x-exp(x) -0.51.52.50.252.5引言上一頁下一頁主 頁 對(duì)于一般的非線性方程，怎樣得到比較精確的近似解，正是該試驗(yàn)要研究或解決的問題。引言上一頁下一頁主 頁方程的常用求解方法1圖形放大法2. 簡(jiǎn)單迭代法3牛頓迭代法返 回上一頁下一頁主 頁 方程 f(x)=0 1）建立坐標(biāo)系，畫曲線f(x)； 2）觀察曲線f(x)

4、與x軸相交的交點(diǎn)； 3）將其中一個(gè)交點(diǎn)進(jìn)行局部放大； 4）該交點(diǎn)的橫坐標(biāo)值就是方程的根。 1圖形放大法上一頁下一頁主 頁例: 求方程 x5 2x2 + 4 = 0 的一個(gè)根.畫方程曲線圖（畫方程曲線圖（tuxfd.m）x=-6:0.01:6;y=x.5+2.*x.2+4;y1=zeros(1,length(x);plot(x,y,x,y1)或或ezplot(f(x),a,b)-6-4-20246-8000-6000-4000-200002000400060008000由此判斷：方程的一個(gè)根在區(qū)間由此判斷：方程的一個(gè)根在區(qū)間-2-2，22內(nèi)，因此將區(qū)間內(nèi)，因此將區(qū)間-6-6，66縮小至縮小至-2

5、-2，22，再觀察！，再觀察！ 該方程有幾個(gè)根？欲尋找其中一個(gè)實(shí)根，并且達(dá)到一定的精度。1圖形放大法上一頁下一頁主 頁 逐次縮小區(qū)間，觀察一個(gè)根在-1.55-1.5之間。MATLAB(tuxfd1.m)-202-50050-2-1.5-1-20-10010-2-1.5-20-10010-1.6-1.5-1.4-2-1011圖形放大法返 回上一頁下一頁主 頁2. 2. 簡(jiǎn)單迭代算法簡(jiǎn)單迭代算法引例：引例： 3xex = 0 Z 1 1）該方程有多少個(gè)根？如何判）該方程有多少個(gè)根？如何判 斷？斷？2 2）如何進(jìn)行迭代求解？）如何進(jìn)行迭代求解？方程變形方程變形： x = ex/3 0 1/3 1/3

6、 0.4652 0.4652 0.5308上一頁下一頁主 頁序號(hào)0123456左邊00.333 0.465 0.531 0.567 0.588 0.599右邊 0.333 0.465 0.531 0.567 0.588 0.599 0.607序號(hào)78910左邊 0.607 0.612 0.615 0.616右邊 0.612 0.615 0.616 0.617方程3xex=0的迭代求解表： x = ex/32. 2. 簡(jiǎn)單迭代算法簡(jiǎn)單迭代算法上一頁下一頁主 頁00.5100.20.40.60.810.616y=xy =ex/3迭代過程如圖所示2. 2. 簡(jiǎn)單迭代算法簡(jiǎn)單迭代算法上一頁下一頁主 頁

7、 方程： f (x) = 0 經(jīng)過簡(jiǎn)單變形：x = j (x) 或 j (x) = f (x)+xx 被稱為不動(dòng)點(diǎn)。迭代過程如下： xn+1 =j (xn)，n =0，1， x0 定義為迭代初值。迭代算法步驟1.1.表達(dá)式表達(dá)式x = j (x)是否唯一是否唯一？2.2.迭代產(chǎn)生的序列是否一定會(huì)收斂？迭代產(chǎn)生的序列是否一定會(huì)收斂？3. 3. 迭代與初始值迭代與初始值x0是否有關(guān)？是否有關(guān)？Z 上一頁下一頁主 頁解：解： 第一步第一步 構(gòu)造迭代函數(shù)： x=j (x)(111)(1)(132232123xxxxxxxxxxxxjjj例：例：用迭代方法求解方程 x3 x2 x1 0。2. 2. 簡(jiǎn)單

8、迭代算法簡(jiǎn)單迭代算法上一頁下一頁主 頁第二步第二步 迭代設(shè)定初值 x0=1， xn+1 =j (xn)，n =0，1， 用 MATLAB 編程（died2.m)x(1)=1;y(1)=1;z(1)=1; %(初始點(diǎn)初始點(diǎn)）for k=1:20 x(k+1)=x(k)3-x(k)2-1; %j j1 1 (x) y(k+1)=(y(k)2+y(k)+1)(1/3);%j j2 2 (y) z(k+1)=1+1/z(k)+1/z(k)2; %j j3 3(z)endx，y，z2. 2. 簡(jiǎn)單迭代算法簡(jiǎn)單迭代算法上一頁下一頁主 頁序號(hào)j3(x)序號(hào)j2(x)j3(x)1j2(x)1.44223.00

9、0081.81751.813621.65371.444491.83851.855431.75322.1716101.83891.829441.79951.6725111.83911.845451.82091.9554121.83921.835561.83081.7730131.83921.841671.83541.8822j1(x)的迭代是失敗的（迭代不收斂 ）。精確解：x=1.8393計(jì)算結(jié)果計(jì)算結(jié)果2. 2. 簡(jiǎn)單迭代算法簡(jiǎn)單迭代算法上一頁下一頁主 頁 迭代函數(shù)j2(x)和j3(x)的選取是成功的。精確解為 x=1.8393。 并且選取函數(shù)j2(x)、j3(x)其收斂速度不一致，前者的速度

10、快些！結(jié)論1、當(dāng)遇到迭代不收斂時(shí)有什么解決辦法？2、如何提高收斂速度？ Z 對(duì)于給定的方程 f(x) = 0, 有多種方式將它改寫成等價(jià)的形式 x = j(x)。但重要的是如何改寫使得序列如何改寫使得序列收斂？收斂？返 回上一頁下一頁主 頁 若 x= j (x) 迭代不收斂，則不直接使用j(x)迭代，而用由j (x)與x的加權(quán)平均，h(x) =l j (x) +(1l)x 進(jìn)行迭代，其中l(wèi)為參數(shù)。顯然 x = h (x) x = j (x) 加速迭代收斂法xn+1 = h (xn) xn+1 = j (xn) 關(guān)鍵是如何確定函數(shù)關(guān)鍵是如何確定函數(shù)h(x) 中的參數(shù)中的參數(shù)l ？上一頁下一頁主

11、頁理論證明：在滿足|h(x)|tol*x(i)&(in) error(n is full); endend y=x(i);其中輸入初值x0，允許最大迭代次數(shù)n及相對(duì)誤差限tol，輸出根的近似值y。程序中需要調(diào)用函數(shù)M文件fun1.m和dfun1.m。上一頁下一頁主 頁算例 求方程 10, 2, 4 .250)arcsin(2(222vrLvxrLxrxLr的根. 建立兩個(gè)M函數(shù)文件:function y=fun2(x)L=25.4;r=2;v=10;y=L*r2*(pi/2)-asin(x/r)-L*x*sqrt(r2-x2)-v; function y=dfun2(x)L=25.4;

12、r=2;v=10;y=L*(2-r)/sqrt(r2-x2)-L*sqrt(r2-x2);上一頁下一頁主 頁算例 在工作空間中使用如下命令在工作空間中使用如下命令:x0=0.5;tol=0.01;n=10000;y=newton1(x0,n,tol)計(jì)算結(jié)果計(jì)算結(jié)果:y = 1.7166 - 0.0050i思考思考: : 牛頓迭代法與基本迭代加速之間的關(guān)系。牛頓迭代法與基本迭代加速之間的關(guān)系。上一頁下一頁主 頁二次函數(shù)為參數(shù)其中axaxxf),1 ()(考慮a 的取值范圍0, 4。 為能了解迭代數(shù)列xk 收斂或發(fā)散的原因，通過幾何圖形觀察使用迭代格式 ： x0=c xk+1 = f(xk)。進(jìn)

13、行迭代的過程。 函數(shù)的迭代是數(shù)學(xué)研究中的一個(gè)非常重要的思想工具。哪怕是對(duì)一個(gè)相當(dāng)簡(jiǎn)單的函數(shù)進(jìn)行迭代，都可以產(chǎn)生異常復(fù)雜的行為，并由此衍生了一些嶄新的學(xué)科分支，如分形與混沌。例如迭代過程的圖解上一頁下一頁主 頁4 , 0),1 (1axaxxkkk注意：不同的參數(shù)a，產(chǎn)生的序列xn可能收斂或可能發(fā)散。觀察a取什么值時(shí)收斂（或發(fā)散）？三種圖形顯示方式：1）線性連接圖線性連接圖 用點(diǎn)(n, xn)描述迭代點(diǎn)，并用直 線連接這些點(diǎn)所形成的折線圖（橫坐標(biāo)表示n，縱坐標(biāo)表示 xn）。我們可以用圖形的方式來反映迭代的過程。如：迭代過程的圖解上一頁下一頁主 頁震蕩發(fā)散情形收斂情形0102000.10.20.3

14、0.4a=0.5,x0=0.5010200.20.40.60.81a=2.5,x0=0.5010200.20.40.60.81a=3.1,x0=0.5010200.20.40.60.81a=3.5,x0=0.5線性連接圖上一頁下一頁主 頁a=0.5;x1=;x1(1)=0.5;for i=2:20 x1(i)=a*x1(i-1)*(1-x1(i-1);endn=1:20;subplot(2,2,1),plot(n,x1),title(a=0.5,x0=0.5)subplot(2,2,2),plot(n,x2),title(a=1.5,x0=0.5)subplot(2,2,3),plot(n,x

15、3),title(a=2.5,x0=0.5)subplot(2,2,4),plot(n,x4),title(a=3.5,x0=0.5)MATLAB 程序程序died4.m上一頁下一頁主 頁 在直角坐標(biāo)系中，首先畫出直線y= x和曲線y= f(x)，其中f(x)為迭代曲線。 從直線y= x上點(diǎn)An(xn,xn)到曲線y= f(x)上點(diǎn)Bn(xn,xn+1)；（垂直） 從Bn(xn,xn+1)到An(xn+1,xn+1); (水平） 重復(fù)。a=2.9, x0=0.2 -0.200.20.40.60.811.2-0.200.20.40.60.81AnBnAny = xy = f(x)蛛網(wǎng)圖上一頁下一

16、頁主 頁發(fā)散情形收斂情形00.51-0.200.20.40.60.811.2xxa=2.5, x0=0.1 00.51-0.200.20.40.60.811.2xxa=1.5, x0=0.1 00.51-0.200.20.40.60.811.2xxa=3.4, x0=0.2 00.51-0.200.20.40.60.811.2xxa=3.7, x0=0.2, n=80 蛛網(wǎng)圖上一頁下一頁主 頁x1=;a=2.9;x=-0.2:0.01:1.2;y=a.*x.*(1-x);plot(x,y),hold on % 畫二次函數(shù)曲線ezplot(x,-0.2,1.2)% 畫直線x1(1)=0.2; %

17、 初始點(diǎn)for i=2:50 x1(i)=a*x1(i-1)*(1-x1(i-1); plot(x1(i-1),x1(i-1),x1(i-1),x1(i); plot(x1(i-1),x1(i),x1(i),x1(i);End % 畫折線x1;died8.mMATLAB 程序程序上一頁下一頁主 頁該圖形特點(diǎn)：觀察當(dāng)參數(shù)a連續(xù)變化時(shí)，函數(shù)f(x)=ax(1-x)的收斂與發(fā)散情況。a=2.9a=3.4混沌圖上一頁下一頁主 頁 在連接圖和蛛網(wǎng)圖中，可以觀察對(duì)不同參數(shù)參數(shù)值，函數(shù)迭代序列的斂散性情況。但參數(shù)僅取幾個(gè)離散值，導(dǎo)致序列收斂與發(fā)散的參數(shù)a的分界值究竟在什么位置？能否對(duì)參數(shù)連續(xù)取值觀察序列收斂

18、性的變化？axn a1a2混沌圖的構(gòu)圖原理上一頁下一頁主 頁 具體做法是：將區(qū)間（0，4以某個(gè)步長(zhǎng)a離散化。對(duì)每個(gè)離散的a值做迭代，忽略前50個(gè)迭代值x1,x2,x50，而把點(diǎn)(a, x51), (a, x52), , (a, x100)顯示在坐標(biāo)平面上，最后形成的圖形稱為Feigenbaum圖。該圖主要目的是：對(duì)不同的參數(shù)a系統(tǒng)地觀察迭代的行為?；煦鐖D的構(gòu)圖原理上一頁下一頁主 頁編寫一個(gè)對(duì)含參變量函數(shù) f(x,a) 進(jìn)行迭代可調(diào)用程序。function root=iter(x0,a)x=;x(1)=x0; for i=2:100 x(i)=a*x(i-1)*(1-x(i-1);endroot

19、=x; % 產(chǎn)生100個(gè)迭代序列組成的數(shù)組MATLAB 程序程序iter.m上一頁下一頁主 頁clf; x=2;hold on;for a=2:0.01:4 root=iter(x,a); plot(a.*ones(size(root(51:100),root(51:100),.)endxlabel(parameter a);ylabel(迭代序列（迭代序列（51-10051-100）); died9.mMATLAB 程序程序上一頁下一頁主 頁 從極限分支點(diǎn)之后，F(xiàn)eigenbaum圖顯得很雜亂，似乎沒有任何規(guī)律。實(shí)際上，對(duì)任何初始值做迭代都會(huì)得到同樣的結(jié)果。這就是所謂的混沌現(xiàn)象。 我們以a=

20、4為例來說明迭代序列對(duì)初值的敏感性。 觀察：無論兩個(gè)初值如何接近，在迭代過程中，與之產(chǎn)生的對(duì)應(yīng)序列xn、yn將漸漸分開，稱迭代對(duì)初值是敏感的。這是任何混沌系統(tǒng)都具有的特性之一。這種特性使得混沌系統(tǒng)會(huì)產(chǎn)生似乎是隨機(jī)的沒有規(guī)律的現(xiàn)象?；煦鐖D的特點(diǎn)上一頁下一頁主 頁迭代格式： xn+1=4xn(1-xn), n = 0,1,1、取初始點(diǎn)x0=0.2, 記產(chǎn)生的序列為xn；2、取初始點(diǎn)x0=0.201, 記產(chǎn)生的序列為yn；3、作線性連接圖。即 plot(n, xn-yn)；01020304050-1-0.8-0.6-0.4-0.200.20.4a=4,x0=0.5,x0=0.510102030405

21、0-1-0.500.51a=4,x0=0.2,x0=0.201敏感性分析上一頁下一頁主 頁 另外，混沌不是隨機(jī)的，由圖形可以發(fā)現(xiàn)，它們具有自相似現(xiàn)象。呈周期變化的形態(tài)。01020304050-1-0.8-0.6-0.4-0.200.20.4a=4,x0=0.5,x0=0.5101020304050-1-0.500.51a=4,x0=0.2,x0=0.201敏感性分析上一頁下一頁主 頁 通過實(shí)驗(yàn)我們發(fā)現(xiàn)函數(shù)f(x)=ax(1-x)對(duì)參數(shù)a是敏感的。當(dāng)參數(shù)取值a = 2.9時(shí)，迭代是收斂的，且收斂于0.655；而對(duì)a = 3.4時(shí)，迭代是發(fā)散的，且在兩個(gè)值之間有規(guī)律地來回取值，因此稱為2-周期的，

22、當(dāng)a = 3.55左右，其迭代是4-周期的；然后是8-周期，16-周期，等等。這就是所謂的分歧現(xiàn)象。隨著參數(shù)a取值的增加，分歧頻率加快。費(fèi)根鮑通過計(jì)算機(jī)進(jìn)行數(shù)值計(jì)算驚奇地發(fā)現(xiàn)：分歧頻率是成比例增大，其比例常數(shù)是4.6692016。 小 結(jié)返 回上一頁下一頁主 頁迭代迭代分形分形（1）（2）（3）（4）上一頁下一頁主 頁迭代迭代分形分形雪花曲線生成元00.10.20.30.40.50.60.70.80.9100.050.10.150.20.250.30.3500.10.20.30.40.50.60.70.80.9100.050.10.150.20.250.30.35Koch.m上一頁下一頁主 頁

23、迭代迭代分形分形樹木花草生成元上一頁下一頁主 頁迭代迭代分形分形Minkowski曲線生成元上一頁下一頁主 頁迭代迭代分形分形Sierpinski三角形生成元012345678012345670123456780123456701234567801234567Serpinski.m01234567801234567上一頁下一頁主 頁【問題背景】【問題背景】 一段時(shí)間一段時(shí)間, 美國(guó)原子能委員會(huì)是按以下方式處理濃縮放射美國(guó)原子能委員會(huì)是按以下方式處理濃縮放射性廢物的性廢物的. 他們將廢物裝入密封性能很好的圓桶中他們將廢物裝入密封性能很好的圓桶中, 然后扔到水然后扔到水深深300英尺的海里英尺的海

24、里. 這種做法是否會(huì)造成放射性污染這種做法是否會(huì)造成放射性污染, 很自然地很自然地引起了生態(tài)學(xué)家及社會(huì)各界的關(guān)注引起了生態(tài)學(xué)家及社會(huì)各界的關(guān)注. 原子能委員會(huì)一再保證原子能委員會(huì)一再保證, 圓圓桶非常堅(jiān)固桶非常堅(jiān)固, 決不會(huì)破漏決不會(huì)破漏, 這種做法是絕對(duì)安全的這種做法是絕對(duì)安全的. 然而一些工然而一些工程師們卻對(duì)此表示懷疑程師們卻對(duì)此表示懷疑, 他們認(rèn)為圓桶在海底相撞時(shí)有可能發(fā)他們認(rèn)為圓桶在海底相撞時(shí)有可能發(fā)生破裂生破裂. 由此雙方展開了一場(chǎng)筆墨官司由此雙方展開了一場(chǎng)筆墨官司. 究竟誰的意見正確呢究竟誰的意見正確呢? 只能讓事實(shí)說話了只能讓事實(shí)說話了!范例: 放射性廢物的處理問題上一頁下一頁

25、主 頁 問題的關(guān)鍵在于圓桶到底能承受多大速度問題的關(guān)鍵在于圓桶到底能承受多大速度的碰撞的碰撞? 圓桶和海底碰撞時(shí)的速度有多大圓桶和海底碰撞時(shí)的速度有多大? 工程師們進(jìn)行了大量破壞性的實(shí)驗(yàn)工程師們進(jìn)行了大量破壞性的實(shí)驗(yàn), 發(fā)現(xiàn)發(fā)現(xiàn)圓桶在直線速度為圓桶在直線速度為40 ft/s 的沖撞下會(huì)發(fā)生破裂的沖撞下會(huì)發(fā)生破裂, 剩下的問題就是計(jì)算圓桶沉入剩下的問題就是計(jì)算圓桶沉入300 ft 深的海底深的海底時(shí)時(shí), 其末速度究竟有多大其末速度究竟有多大? 問題分析問題分析范例上一頁下一頁主 頁 1. 使用使用55加侖的圓桶加侖的圓桶; ( 1加侖加侖 = 3.7854升升 ) 2. 裝滿放射性廢物時(shí)的圓桶重量為裝滿放射性廢物時(shí)的圓桶重量為 W = 527.436磅磅 (1 磅磅 = 0.4526公斤公斤 ) 3. 在海水中圓桶受到的浮力在海水