數(shù)學(xué)實(shí)驗(yàn)第三次講稿

1、2022-4-131 1) 掌握方程求解的三種解法：解析 法、數(shù)值解法以及圖形表示解的方法; 2) 學(xué)會(huì)使用MATLAB軟件求解析解、數(shù)值解 和圖形解； 3） 通過(guò)范例學(xué)習(xí)怎樣建立方程模型和 分析問(wèn)題的思想。實(shí)驗(yàn)?zāi)康?022-4-132一，引例 例子【問(wèn)題背景問(wèn)題背景】一段時(shí)間一段時(shí)間, 美國(guó)原子能委員會(huì)是美國(guó)原子能委員會(huì)是按以下方式處理濃縮放射性廢物的按以下方式處理濃縮放射性廢物的. 他們將廢物裝入他們將廢物裝入密封性能很好的圓桶中密封性能很好的圓桶中, 然后扔到水深然后扔到水深300英尺的海英尺的海里里. 這種做法是否會(huì)造成放射性污染這種做法是否會(huì)造成放射性污染, 很自然地引起很自然地引起

2、了生態(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)為他們認(rèn)為圓桶在海底相撞時(shí)有可能發(fā)生破裂圓桶在海底相撞時(shí)有可能發(fā)生破裂. 由此雙方展開了由此雙方展開了一場(chǎng)筆墨官司一場(chǎng)筆墨官司. 究竟誰(shuí)的意見正確呢究竟誰(shuí)的意見正確呢? 只能讓事實(shí)只能讓事實(shí)說(shuō)話了說(shuō)話了!2022-4-133數(shù)學(xué)建模的一般步驟模型準(zhǔn)備模型準(zhǔn)備模型假設(shè)模型假設(shè)模型構(gòu)成模型構(gòu)成模型檢驗(yàn)?zāi)Ｐ蜋z驗(yàn)?zāi)Ｐ头治瞿Ｐ头治瞿Ｐ颓蠼饽Ｐ?/p>

3、求解模型應(yīng)用模型應(yīng)用2022-4-134二，方程求解 1，解析方法 2，圖形放大法 3，迭代方法 4，區(qū)間方法2022-4-1351，方程求解之解析方法 主要針對(duì)一些比較簡(jiǎn)單的方程以及方程組，比如多項(xiàng)式方程等。同學(xué)們以前對(duì)方程的求解也是針對(duì)這樣一些方程進(jìn)行的。該方法的優(yōu)點(diǎn)是可以利用紙筆得到簡(jiǎn)單有效并且精確的解；缺點(diǎn)是可以求解的方程數(shù)量太少。Matlab和Maple提供了求方程解析解的函數(shù)，可以說(shuō)對(duì)數(shù)學(xué)演算提供了不少方便。2022-4-1362，方程求解之圖形放大法 圖形最大的有點(diǎn)就是直觀，試想如果我們有了函數(shù)準(zhǔn)確的圖形，那么曲線和x軸的交點(diǎn)就是我們要求的方程的解。因此我們可以利用圖形工具得到方

4、程的解。當(dāng)然，計(jì)算機(jī)上的圖形不可能等同于函數(shù)的真實(shí)圖形，因?yàn)橛?jì)算機(jī)上的圖形是曲線上部分點(diǎn)的軌跡而不是全部，因此通過(guò)圖形不可能得到方程的精確解，甚至它只是一個(gè)比較粗略的解，當(dāng)然，通過(guò)對(duì)圖形的放大可以得到更精確一些的解。同時(shí)，這種方法也不適應(yīng)大量的數(shù)據(jù)處理。2022-4-137 方程 f(x)=0 1）建立坐標(biāo)系，畫曲線f(x)； 2）觀察曲線f(x)與x軸相交的交點(diǎn)； 3）將其中一個(gè)交點(diǎn)進(jìn)行局部放大； 4）該交點(diǎn)的橫坐標(biāo)值就是方程的根。 2022-4-138例: 求方程 x5 +2x2 + 4 = 0 的一個(gè)根.畫方程曲線圖（畫方程曲線圖（tuxfd.m）x=-3:0.01:3;y=x.5+2*

5、x.2+4;y1=0*x;plot(x,y,x,y1)-6-4-20246-8000-6000-4000-200002000400060008000由此判斷：方程的一個(gè)根在區(qū)間由此判斷：方程的一個(gè)根在區(qū)間-2-2，22內(nèi)，因此將區(qū)間內(nèi)，因此將區(qū)間-3-3，33縮小至縮小至-2-2，22，再觀察！，再觀察！ 該方程有幾個(gè)根？欲尋找其中一個(gè)實(shí)根，并且達(dá)到一定的精度。2022-4-139 逐次縮逐次縮小區(qū)間，觀小區(qū)間，觀察一個(gè)根在察一個(gè)根在-1.55-1.5之之間。間。-202-50050-2-1.5-1-20-10010-2-1.5-20-10010-1.6-1.5-1.4-2-1012022-4

6、-13103,方程求解之迭代法 迭代法的理論以及方法的出現(xiàn)，對(duì)方程求解有著里程碑式的意義。其基本思想如下: 需要求解的方程： f (x) = 0 （1） 經(jīng)過(guò)某種變形得：x = j (x) （2） 從而求解方程（1）轉(zhuǎn)化成為求解（2）得不動(dòng)點(diǎn)。（滿足條件x*=j(x*)的點(diǎn)x*稱為不動(dòng)點(diǎn)） 為得到方程的不動(dòng)點(diǎn)，可以構(gòu)造迭代過(guò)程如下： xn+1 =j (xn)，n =0，1， x0 定義為迭代初值。2022-4-1311解：解： 第一步第一步 構(gòu)造迭代函數(shù)： x=j (x)x(x1x11x)x(1xxx)x(1xxx32232123j+j+j例：例：用迭代方法求解方程 x3 x2 x1 0。20

7、22-4-1312第二步第二步 迭代設(shè)定初值 x0=1， xn+1 =j (xn)，n =0，1， 用 MATLAB 編程（died2.m文件)x=1;y=1;z=1;(初始點(diǎn)）for k=1:20 x=x3-x2-1; % j1 (x) y=(y2+y+1)(1/3); % j2 (y) z=1+1/z+1/z2; % j3 (z)endX，y，z2022-4-1313序號(hào)j3(x)序號(hào)j2(x)j3(x)1j2(x)1.44223.000081.81751.813621.65371.444491.83851.855431.75322.1716101.83891.829441.79951.6

8、725111.83911.845451.82091.9554121.83921.835561.83081.7730131.83921.841671.83541.8822j1(x)的迭代是失敗的（迭代不收斂 ）。精確解：x=1.8393計(jì)算結(jié)果計(jì)算結(jié)果2022-4-1314 迭代函數(shù)j2(x)和j3(x)的選取是成功的。精確解為 x=1.8393。 并且選取函數(shù)j2(x)、j3(x)其收斂速度不一致，前者的速度快些！結(jié)論1、當(dāng)遇到迭代不收斂時(shí)有什么解決辦法？2、如何提高收斂速度？ 對(duì)于給定的方程 f(x) = 0, 有多種方式將它改寫成等價(jià)的形式 x = j(x)。但重要的是如何改寫使得序列如何

9、改寫使得序列收斂？收斂？2022-4-1315 當(dāng)今最流行的迭代法是牛頓法以及由此改進(jìn)的一些方法，比如擬牛頓法等。其基本的思想就是構(gòu)造迭代格式是利用函數(shù)的導(dǎo)數(shù)，這類方法有收斂速度快，穩(wěn)定性好等特點(diǎn)。對(duì)低維和高維情況都適合，也是當(dāng)今一些軟件均采用的方法。當(dāng)然，因?yàn)樾枰瘮?shù)的導(dǎo)數(shù)信息，所以自然對(duì)不可微的問(wèn)題受到制約。該方法的迭代格式為： xk+1=xk-f(xk)/f(xk)2022-4-1316例:利用牛頓法求方程x3 -x2 -x-1 = 0的根.方法:第一步,給出函數(shù)的導(dǎo)函數(shù)3x2-2x-1; 第二步,給出函數(shù)的迭代格式: xk+1=xk-f(xk)/f(xk); 設(shè)置一定的精度要求,達(dá)到即

10、終止.定義函數(shù)m文件: function ff=mynewton(x) ff=(x3 x2 -x-1 )/(3*x2-2*x-1);定義命令m文件: x0=1; x1=x0-mynewton (x0); while abs(x1-x0)0.0001 x0=x1; x1=x0-mynewton (x0); end2022-4-1317三，解方程函數(shù)格式及例子 Matlab對(duì)方程的求解提供了以下的一些函數(shù)： （1）多項(xiàng)式求根； （2） 線性方程組求解； （3）一般的非線性方程（組）求解：2022-4-1318-1.3-1.2-1.1-1-1.5-1-0.500.51x0輸出: -1.2131 -0

11、.9017 + 0.5753i -0.9017 - 0.5753i -0.2694 + 0.9406i -0.2694 - 0.9406i 0.4168 + 0.8419i 0.4168 - 0.8419i 0.8608 + 0.3344i 0.8608 - 0.3344i例：求解多項(xiàng)式方程 x9+x8+1=0輸入: p=1,1,0,0,0,0,0,0,0,1; roots(p)1，roots()語(yǔ)句的用法-10-50510-2-1.5-1-0.500.511.52x 108x02022-4-1319線性方程組：AX = bAX = b 其中A是mn階矩陣，b是m維向量。 x=A b or x

12、=inv(A)*b特點(diǎn)：只能求出一個(gè)特解。2022-4-13202，線性方程組solution to the following linear system of equations: You can formulate and solve the problem as A = 3 11 -2; 1 1 -2; 1 -1 1; b = 7; 4; 19; x = Ab x = 13.2188 -2.3438 3.43752022-4-13213146,987654321Ab例: AX = b, 解: 輸入:A=1 2 3; 4 5 6; 7 8 9; b=6; 14; -3; x1=Ab, x

13、2=inv(A)*b輸出:警告: 系統(tǒng)的秩不足. 解不唯一.1、題中rank(A)=rank(A|b)=23, 該方程組有無(wú)窮解。2、輸出結(jié)果是否一致？3、如何求方程組的全部解？思考返 回2022-4-13223.1，函數(shù)fzero的用法格式： fzero(函數(shù)名,初值或區(qū)間) Example 1. Calculate by finding the zero of the sine function near 3. x = fzero(sin,3) x = 3.1416 Example 2. To find the zero of cosine between 1 and 2 x = fzer

14、o(cos,1 2) x = 1.57082022-4-13231、方程(組)， f1(x) = 0,fn(x) = 0, x = (x1,xn)2、方程(組)， f1(x) = 0,fn(x) = 0, x = (x1,xn)fun.m function f = fun(x) f(1)= f1(x) ; f(n)= fn(x) 初值1）可以省略。2）options=1,表示輸出中間結(jié)果。solve(f1(x),f2(x),fn(x) ) X = fsolve (fun, X0, options) MATLAB軟件直接求解法2022-4-1324輸出: 1/2/a*(-b+(b2-4*a*c)

15、(1/2) 1/2/a*(-b-(b2-4*a*c)(1/2)單變量方程單變量方程0)(xf3.2 solve()語(yǔ)句的用法例1： 求解方程 ax2+bx+c = 0輸入: x= solve(a*x2+b*x+c)或 solve(a*x2+b*x+c=0)1）符號(hào)解2022-4-1325例2： 解方程： x3-2x2=x-1解： s=solve(x3-2*x2=x-1) double(s)2）數(shù)字解該方程是否有實(shí)根？vpa(s,10) 2022-4-1326 例3 求解方程： tan(x)-sin(x)=0 3）無(wú)窮解輸入: solve(tan(x)-sin(x)=0)輸出：0 (不能給出全部

16、解） (tx1.m)2022-4-1327輸入： x,y=solve(x2*y2, x-(y/2)-b)輸出：x = 0 ， y = -2*b 0， -2*b （符號(hào)解） b， 0 b， 0 v=x,y 多變量方程組多變量方程組0)(,0)(1xfxfm例4byxyx2022byxyx2122或 b=2solve()語(yǔ)句的用法2022-4-1328+05012307lnsin32zyxzxzyxy例6：求解方程組解 輸入：syms x,y,z;x,y,z=solve(sin(x)+y2+log(z)-7=0, 3*x+2y-z3+1=0,x+y+z-5=0,x,y,z)x =0.5990537

17、5664056731520568183824539y =2.3959314023778168490940003756591z =2.0050148409816158357003177860955輸出：3.3 fsolve()語(yǔ)句的用法2022-4-1329解：1）建立方程組的M-函數(shù)文件(nxxf.m) function eq=nxxf(x) eq(1)=sin(x(1)+x(2)2+log(x(3)-7; eq(2)=3*x(1)+2x(2)-x(3)3+1; eq(3)=x(1)+x(2)+x(3)-5;2) 運(yùn)行程序(test4.m) y=fsolve(nxxf,1,1,1)3）運(yùn)行結(jié)果

18、：Optimization Terminated Successfully y= 0.5990 2.3959 2.0050fsolve()語(yǔ)句的用法2022-4-1330fsolve()函數(shù)的第三個(gè)輸入是函數(shù)的第三個(gè)輸入是options,它是一個(gè)結(jié)構(gòu)型數(shù)據(jù)它是一個(gè)結(jié)構(gòu)型數(shù)據(jù),可以可以通過(guò)函數(shù)通過(guò)函數(shù)optimset ()進(jìn)行設(shè)定進(jìn)行設(shè)定.當(dāng)不進(jìn)行設(shè)定時(shí)采用缺省設(shè)置當(dāng)不進(jìn)行設(shè)定時(shí)采用缺省設(shè)置. fsolve函數(shù)還可以有后面的參數(shù)設(shè)定函數(shù)還可以有后面的參數(shù)設(shè)定,這一功能在有的時(shí)候非常這一功能在有的時(shí)候非常有用有用.fsolve()語(yǔ)句的用法例如,求解函數(shù)sin(ax)-x=0的最小正解.分析: 已

19、知的知識(shí)我們知道,該方程沒有解析形式的解,也就是說(shuō)很難得到這個(gè)解和a之間的具體關(guān)系.數(shù)值的方法可以對(duì)a先取確定的值,這樣變成一個(gè)一元方程,容易進(jìn)行求解,當(dāng)我們變化a時(shí),就得到很多這樣的解,通過(guò)這樣的方法,我們可以得到解和a之間的一些大致關(guān)系.2022-4-1331具體方法如下具體方法如下:建立函數(shù)建立函數(shù)m文件文件:function ff=funpara(x,a)ff=sin(a*x)-x;fsolve()語(yǔ)句的用法建立相應(yīng)的命令m文件:B=zeros(100,1); for a=1:100 x0=pi/(2*a)+0.01;B(a)=fsolve(funpara,x0,a);endplot(

20、B)2022-4-13322022-4-1333例：例：某物體的邊緣呈圓形，通過(guò)測(cè)量邊沿上是一個(gè)點(diǎn)的坐標(biāo)，數(shù)據(jù)為（見程序）。使用fsolve計(jì)算物體邊緣的曲線方程。解：假定物體邊緣的曲線方程為（x-a）2+(y-b)2=r2為了得到a，b和r的值。將上述十一組數(shù)據(jù)代人方程，得到（xi-a）2+(yi-b)2=r2（i=1,11）。次數(shù)方程個(gè)數(shù)是11個(gè)，未知量個(gè)數(shù)是3個(gè)。使用fsolve求解，其函數(shù)m文件為：function y=funcircle(x)A=6.7630 5.1313 2.4713 -0.3435 -2.3887 -2.9927 -1.9572 0.3778 3.2455 5.7

21、042 6.9465;B=23.2879 25.6492 26.7268 26.1668 24.1531 21.3470 18.6699 17.0010 16.8883 18.3688 20.9564;for i=1:11 y(i)=(A(i)-x(1)2+(B(i)-x(2)2-x(3)2;end二分法 對(duì)于一個(gè)閉區(qū)間上的連續(xù)函數(shù)，我們有一個(gè)0點(diǎn)存在定理，利用這個(gè)定理，可以對(duì)不可微的函數(shù)求得函數(shù)的0點(diǎn)。 基本思想是通過(guò)判斷函數(shù)在端點(diǎn)處的函數(shù)值異號(hào)可以確定函數(shù)在開區(qū)間上至少有一個(gè)0點(diǎn)，然后通過(guò)縮小區(qū)間得到解的近似。 優(yōu)點(diǎn)是不需要函數(shù)的導(dǎo)數(shù)信息，而且只要有0點(diǎn)就一定可以得到；缺點(diǎn)是相對(duì)牛頓法等速

22、度較慢。2022-4-1334 問(wèn)題的關(guān)鍵在于圓桶到底能承受多大速度問(wè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ā)生破裂, 剩下的問(wèn)題就是計(jì)算圓桶沉入剩下的問(wèn)題就是計(jì)算圓桶沉入300 ft 深的海底深的海底時(shí)時(shí), 其末速度究竟有多大其末速度究竟有多大? 問(wèn)題分析問(wèn)題分析引例的分析和求解2022-4-1335 1. 使用使用55加侖的圓桶加侖的圓桶; ( 1加侖加侖 = 3.7854

23、升升 ) 2. 裝滿放射性廢物時(shí)的圓桶重量為裝滿放射性廢物時(shí)的圓桶重量為 W = 527.436磅磅 (1 磅磅 = 0.4526公斤公斤 ) 3. 在海水中圓桶受到的浮力在海水中圓桶受到的浮力 B = 470.327磅磅 4. 圓桶下沉?xí)r受到海水的阻力圓桶下沉?xí)r受到海水的阻力 D = C v C 為常數(shù)為常數(shù), 經(jīng)測(cè)算得經(jīng)測(cè)算得: C = 0.08. 5. 建立坐標(biāo)系建立坐標(biāo)系, 取垂直向下為坐標(biāo)方向取垂直向下為坐標(biāo)方向 y , 海平面為坐標(biāo)原點(diǎn)海平面為坐標(biāo)原點(diǎn).y0 問(wèn)題假設(shè)問(wèn)題假設(shè)引例2022-4-1336根據(jù)牛頓第二定律根據(jù)牛頓第二定律, , 圓桶下沉?xí)r應(yīng)滿足微分方程圓桶下沉?xí)r應(yīng)滿足微

24、分方程: :)(/(86.713)(lim)1 ()(0)0(:dd,dd22時(shí)當(dāng)速度：容易計(jì)算出圓桶的極限其解：初值條件其中阻力）浮力（重力tsftCBWtveCBWtvvBCvWdtdvmvtyCvDgWmDBWtymtWCg327.470B 建立模型建立模型引例2022-4-1337 為了求出圓桶與海底的碰撞速度為了求出圓桶與海底的碰撞速度v(t), 需要求出需要求出圓桶下沉到海底圓桶下沉到海底300英尺時(shí)的時(shí)間英尺時(shí)的時(shí)間 t, 再計(jì)算再計(jì)算v(t)，要要做到這一點(diǎn)是十分困難的做到這一點(diǎn)是十分困難的. 若將速度若將速度v 看成是海水深看成是海水深度度y 的函數(shù)的函數(shù), 即即vyvtyyvtvtytyvtvdddddddddd)()(22由復(fù)合函數(shù)的求導(dǎo)法知由復(fù)合函數(shù)的求導(dǎo)法知 建立模型建立模型引例2022-4-1338WgyBWCvBWCBWCvyvWgyvCvBWvCvBWyvmvln0)0(,0)0(dddd2積分，得：初值條件：或微分方程變?yōu)椋航柚鷶?shù)值方