數(shù)值計(jì)算方程式求根篇ppt課件

1、教材-應(yīng)用MATLAB於數(shù)值計(jì)算陳鴻智教材-應(yīng)用MATLAB於數(shù)值計(jì)算u數(shù)值計(jì)算是件很好玩的事，可以用很精簡(jiǎn)的程式完成數(shù)值計(jì)算是件很好玩的事，可以用很精簡(jiǎn)的程式完成u 並且可以將它用圖形的方式展示出來。並且可以將它用圖形的方式展示出來。uMatlab可以讓我們作到。可以讓我們作到。教材-應(yīng)用MATLAB於數(shù)值計(jì)算u數(shù)值計(jì)算介紹數(shù)值計(jì)算介紹u試圖告訴你數(shù)值計(jì)算是件很棒試圖告訴你數(shù)值計(jì)算是件很棒的事的事教材-應(yīng)用MATLAB於數(shù)值計(jì)算u工程需求工程需求u 資料探索分析資料探索分析u “ “透過觀察讓科學(xué)大幅進(jìn)展透過觀察讓科學(xué)大幅進(jìn)展”牛頓牛頓u “ “數(shù)學(xué)就是要求嚴(yán)密、美感和結(jié)構(gòu)數(shù)學(xué)就是要求嚴(yán)密、

2、美感和結(jié)構(gòu)”u - - 高斯高斯牛頓牛頓高斯高斯教材-應(yīng)用MATLAB於數(shù)值計(jì)算uGeorge Forsythe 於於1960年代後期在年代後期在 Stanford 大學(xué)開創(chuàng)了數(shù)值方法的課程。大學(xué)開創(chuàng)了數(shù)值方法的課程。uForsythe, Malcolm和和 Moler合寫了一合寫了一本教材。本教材。u那時(shí)的課程都是以那時(shí)的課程都是以Fortran語言編寫語言編寫的程式庫為主。的程式庫為主。u為解決線性代數(shù)課程中矩陣運(yùn)算為解決線性代數(shù)課程中矩陣運(yùn)算問題，問題， Cleve Moler博士於博士於1980年開年開發(fā)出一套簡(jiǎn)單易用的矩陣運(yùn)算工具，發(fā)出一套簡(jiǎn)單易用的矩陣運(yùn)算工具， Matlab於是誕

3、生。於是誕生。Cleve Moler教材-應(yīng)用MATLAB於數(shù)值計(jì)算 由由 Cleve Moler 研發(fā)將研發(fā)將EISPACK & LINPACK包裝成包裝成一種容易操控的介面一種容易操控的介面1957: Moler到到Caltech. 研究數(shù)值分析研究數(shù)值分析1961:到史丹佛大學(xué)講演到史丹佛大學(xué)講演. 與與 G. Forsythe 合作研究合作研究Laplacian eigenvalues.1977: 第第1版版 Matlab; 2,000條條Fortran語法程式碼語法程式碼80 個(gè)函數(shù)個(gè)函數(shù) (現(xiàn)已超過現(xiàn)已超過 8000個(gè)函數(shù)個(gè)函數(shù) )1979: 在史丹佛與在史丹佛與 Jack

4、 Little 相遇，開始將源碼轉(zhuǎn)換成相遇，開始將源碼轉(zhuǎn)換成C1984: Mathworks 成立成立教材-應(yīng)用MATLAB於數(shù)值計(jì)算uMatlab, Simulink, Stateflowu Matlab 2009bu 各種領(lǐng)域的運(yùn)用各種領(lǐng)域的運(yùn)用u航太航太,國防國防, 電腦電腦,通信通信, 測(cè)繪測(cè)繪u Mathworks 公司公司u3,500 以上的大學(xué)使用以上的大學(xué)使用,500,000以上的用戶以上的用戶,100以上國家使用以上國家使用u2019年收益年收益: 3億美元以上億美元以上.u2019員工員工: 2,000人以上人以上u 價(jià)錢價(jià)錢:u 2,000$ (商務(wù)版價(jià)商務(wù)版價(jià)),u 1

5、00$ (學(xué)生版價(jià)學(xué)生版價(jià))教材-應(yīng)用MATLAB於數(shù)值計(jì)算- 505- 4- 2024- 1- 0 . 500 . 51u研發(fā)單位研發(fā)單位u老師老師u實(shí)驗(yàn)室軟體設(shè)備之一實(shí)驗(yàn)室軟體設(shè)備之一u Matlab 著重在演算法而不是程式語言探討著重在演算法而不是程式語言探討u學(xué)生學(xué)生u批次處理檔案批次處理檔案u不必再需要另外學(xué)其他的語言不必再需要另外學(xué)其他的語言!u給你一個(gè)測(cè)試概念的好環(huán)境給你一個(gè)測(cè)試概念的好環(huán)境u把概念很快轉(zhuǎn)換成程式碼把概念很快轉(zhuǎn)換成程式碼, 然後轉(zhuǎn)成然後轉(zhuǎn)成 C或或Java 等等u易懂的視覺圖易懂的視覺圖教材-應(yīng)用MATLAB於數(shù)值計(jì)算u微積分微積分u常微分方程常微分方程u矩陣矩陣

6、u計(jì)算機(jī)程式編寫經(jīng)驗(yàn)計(jì)算機(jī)程式編寫經(jīng)驗(yàn)教材-應(yīng)用MATLAB於數(shù)值計(jì)算u是求實(shí)數(shù)自變量、實(shí)值函數(shù)是求實(shí)數(shù)自變量、實(shí)值函數(shù)f(x)=0f(x)=0的解之較慢的解之較慢但很考靠的算法但很考靠的算法uf(x)f(x)是個(gè)連續(xù)函數(shù)，且在是個(gè)連續(xù)函數(shù)，且在a,ba,b區(qū)間內(nèi)會(huì)改變區(qū)間內(nèi)會(huì)改變正負(fù)號(hào)，則可以找到一個(gè)正負(fù)號(hào)，則可以找到一個(gè)x x* *，使得，使得f(xf(x* *)=0)=0uMATLABMATLABuk = 0;k = 0;u while abs(b-a) eps while abs(b-a) eps* *abs(b)abs(b)u x = (a + b)/2; x = (a + b)/2

7、;u if sign(f(x) = sign(f(b) if sign(f(x) = sign(f(b)u b = x; b = x;u else elseu a = x; a = x;u end endu k = k + 1; k = k + 1;u end end教材-應(yīng)用MATLAB於數(shù)值計(jì)算u公式xn+1= xn-f(xn)/f(xn)uMATLABu k = 0;u while abs(x - xprev) eps*abs(x)u xprev = x;u x = x - f(x)/fprime(x)u k = k + 1;u endu以計(jì)算為 例u while abs(x - xpr

8、ev) eps*abs(x)u xprev = x;u x = 0.5*(x + M/x);u endu缺點(diǎn)u當(dāng)f(x)不連續(xù)、無有界的一階、二皆導(dǎo)數(shù)，或者初始值不接近準(zhǔn)確解時(shí)就無法使用本方法無限循環(huán)的牛頓法需要推導(dǎo)公式教材-應(yīng)用MATLAB於數(shù)值計(jì)算u用兩次迭代解構(gòu)造出的有限差分近似，替代牛頓法中的求導(dǎo)u公式uMATLABu while abs(b - a) eps*abs(b)u c = a;u a=b;u b = b + (b-c)/f(c)/f(b)-1);u k = k + 1;u endu優(yōu)點(diǎn)u不需要顯式計(jì)算f(x) ，而有類似收斂的性質(zhì)。en+1=O(enen-1)，是超線性收斂

9、，事實(shí)上en+1=O(en)u=(1+5)/2111()()(),nnnnnnnnnf xf xf xsxxxxs教材-應(yīng)用MATLAB於數(shù)值計(jì)算u既然割線法使用前兩個(gè)近似解得到下一個(gè)解，那麼何既然割線法使用前兩個(gè)近似解得到下一個(gè)解，那麼何不利用前三個(gè)近似解呢不利用前三個(gè)近似解呢? ?u作法：作法：u用三個(gè)點(diǎn)用三個(gè)點(diǎn)a,b,ca,b,c及其函數(shù)值及其函數(shù)值f(a), f(b), f(c)f(a), f(b), f(c)，組一，組一個(gè)個(gè)2 2次曲線次曲線u對(duì)數(shù)據(jù)進(jìn)行插值，然後另其與對(duì)數(shù)據(jù)進(jìn)行插值，然後另其與x x軸的交點(diǎn)為下一個(gè)迭代軸的交點(diǎn)為下一個(gè)迭代點(diǎn)點(diǎn)uMATLABMATLABu k=0;

10、k=0;u while abs(c - b) eps while abs(c - b) eps* *abs(c)abs(c)u x = ployinterp(f(a),f(b),fc,a,b,c,0) x = ployinterp(f(a),f(b),fc,a,b,c,0)u a=b; a=b;u b = c b = cu c=x; c=x;u k = k + 1; k = k + 1;u end endu缺點(diǎn)缺點(diǎn)u速度不穩(wěn)定，需要好好馴服速度不穩(wěn)定，需要好好馴服( (人工多次反復(fù)測(cè)試人工多次反復(fù)測(cè)試) )教材-應(yīng)用MATLAB於數(shù)值計(jì)算u結(jié)合二分法的可靠性與割線法及逆二次差值法結(jié)合二分法的可

11、靠性與割線法及逆二次差值法(IQI)(IQI)算法的收斂速度算法的收斂速度u演算法：演算法：u選取初始值選取初始值a a和和b b，使得，使得f(a)f(a)和和f(b)f(b)的正負(fù)號(hào)正好相反的正負(fù)號(hào)正好相反u使用一步割線法，得到使用一步割線法，得到a a和和b b之間的一個(gè)值之間的一個(gè)值u重復(fù)下面的步驟，直到重復(fù)下面的步驟，直到|b -a|b|b -a|function b = fzerotx(F,ab,varargin) if ischar(F) & exist(F)=2 F = inline(F); elseif isa(F,sym) F = inline(char(F); e

12、nd a = ab(1); b = ab(2); fa = feval(F,a,varargin:); fb = feval(F,b,varargin:); if sign(fa) = sign(fb) error(Function must change sign on the interval) end c = a; fc = fa; d = b - c; e = d; 教材-應(yīng)用MATLAB於數(shù)值計(jì)算 while fb = 0 if sign(fa) = sign(fb) a = c; fa = fc; d = b - c; e = d; end if abs(fa) c = b; b

13、= a; a = c; fc = fb; fb = fa; fa = fc; end m = 0.5*(a - b); tol = 2.0*eps*max(abs(b),1.0); if (abs(m) break end if (abs(e) tol) | (abs(fc) d = m; e = m; else s = fb/fc; 教材-應(yīng)用MATLAB於數(shù)值計(jì)算 if (a = c) p = 2.0*m*s; q = 1.0 - s; else q = fc/fa; r = fb/fa; p = s*(2.0*m*q*(q - r) - (b - c)*(r - 1.0); q = (q

14、 - 1.0)*(r - 1.0)*(s - 1.0); end; if p 0, q = -q; else p = -p; end; if (2.0*p 3.0*m*q - abs(tol*q) & (p e = d; d = p/q; else d = m; e = m; end; end 教材-應(yīng)用MATLAB於數(shù)值計(jì)算 c = b; fc = fb; if abs(d) tol b = b + d; else b = b - sign(b-a)*tol; end fb = feval(F,b,varargin:); end教材-應(yīng)用MATLAB於數(shù)值計(jì)算uMATLABufor

15、n = 1:10u z(n) = fzerotx(besselj(0,x),(n-1) n*pi);uendu x = 0:pi/50:10*pi;u y = besselj(0,x);u plot(z,zeros(1,10),o,x,y,-)u line(0 10*pi,0 0,color,black)u axis(0 10*pi -0.5 1.0)輸入:基本函數(shù)、搜索區(qū)間教材-應(yīng)用MATLAB於數(shù)值計(jì)算二分點(diǎn)IQI點(diǎn)割線點(diǎn)教材-應(yīng)用MATLAB於數(shù)值計(jì)算1.給定一個(gè)函數(shù)F(x)和值，求使得F()=;2.給定對(duì)未知函數(shù)F(x)採樣得到的一些數(shù)據(jù)點(diǎn)(xk,yk)和值，求使得F()=;解：1.

16、f(x)=F(x)-=0，可以求得x=2. f(x)=P(x)-， P(x)是某個(gè)插值函數(shù)(可以用pchiptx(xk,yk,x)或splinetc(xk,yk,x)得到P(x) ，其缺點(diǎn)為：計(jì)算代價(jià)大，需要反復(fù)計(jì)算插值。使用反復(fù)插值法，用的是pchip和spline算法，將xk和yk角色互調(diào)yk需具單調(diào)性，或者至少yk的某個(gè)包含目標(biāo)值的子集單調(diào)。生出另一個(gè)分段多項(xiàng)式Q(yk)=xk。y=，=Q(y)。教材-應(yīng)用MATLAB於數(shù)值計(jì)算uf(x),和它所在的區(qū)間a,b ，求x使得f(x)為極小(或極大)u作法:u黃金分割搜索+拋物線法u區(qū)間步長(zhǎng)會(huì)以-10.618的比例u縮小，經(jīng)過-52/log2

17、(-1)75步u後，區(qū)間步長(zhǎng)大致減小為原始長(zhǎng)度的eps倍。u經(jīng)過開始的一些步後，通常有足夠的歷史訊息，給出區(qū)間內(nèi)3個(gè)不同的點(diǎn)，以及對(duì)應(yīng)的函數(shù)值。u如果這3個(gè)點(diǎn)產(chǎn)生的拋物線的極小(大)值落在區(qū)間內(nèi)，那麼下一點(diǎn)通常選擇這個(gè)極小點(diǎn)，而不是黃金分割點(diǎn)。2131011(35) / 20.382 教材-應(yīng)用MATLAB於數(shù)值計(jì)算u求 的極大值uF = inline(-humps(x);ufmintx(F,-1,2,1.e-4)uStep x f(x)uinit: 0.1458980337 -25.2748253202ugold: 0.8541019662 -20.9035150009ugold: -0.2

18、917960675 2.5391843579upara: 0.4492755129 -29.0885282699upara: 0.4333426114 -33.8762343193upara: 0.3033578448 -96.4127439649ugold: 0.2432135488 -71.7375588319upara: 0.3170404333 -93.8108500149upara: 0.2985083078 -96.4666018623upara: 0.3003583547 -96.5014055840upara: 0.3003763623 -96.5014085540upara: 0.3003756221 -96.50140856032211( )(0.3)0.01(0.9)0.04h xxx極大值教材-應(yīng)用MATLAB於數(shù)值計(jì)算uR. P. Brent, Algorithms for Minimization Without Derivatives, Prentice Hall,Englewood Cliffs, 1973.uG. Dahlquest and A. Bjorck, Numerical Methods, Prentice Hall, Englewood Cliffs, 1974.uT. J. Dekker, Finding a Zero