1、教材-應用MATLAB於數(shù)值計算陳鴻智教材-應用MATLAB於數(shù)值計算u數(shù)值計算是件很好玩的事，可以用很精簡的程式完成數(shù)值計算是件很好玩的事，可以用很精簡的程式完成u 並且可以將它用圖形的方式展示出來。並且可以將它用圖形的方式展示出來。uMatlab可以讓我們作到?？梢宰屛覀冏鞯健＝滩?應用MATLAB於數(shù)值計算u數(shù)值計算介紹數(shù)值計算介紹u試圖告訴你數(shù)值計算是件很棒試圖告訴你數(shù)值計算是件很棒的事的事教材-應用MATLAB於數(shù)值計算u工程需求工程需求u 資料探索分析資料探索分析u “ “透過觀察讓科學大幅進展透過觀察讓科學大幅進展”牛頓牛頓u “ “數(shù)學就是要求嚴密、美感和結構數(shù)學就是要求嚴密、

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

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

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

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

6、u計算機程式編寫經(jīng)驗計算機程式編寫經(jīng)驗教材-應用MATLAB於數(shù)值計算u是求實數(shù)自變量、實值函數(shù)是求實數(shù)自變量、實值函數(shù)f(x)=0f(x)=0的解之較慢的解之較慢但很考靠的算法但很考靠的算法uf(x)f(x)是個連續(xù)函數(shù)，且在是個連續(xù)函數(shù)，且在a,ba,b區(qū)間內(nèi)會改變區(qū)間內(nèi)會改變正負號，則可以找到一個正負號，則可以找到一個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教材-應用MATLAB於數(shù)值計算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以計算為 例u while abs(x - xpr

8、ev) eps*abs(x)u xprev = x;u x = 0.5*(x + M/x);u endu缺點u當f(x)不連續(xù)、無有界的一階、二皆導數(shù)，或者初始值不接近準確解時就無法使用本方法無限循環(huán)的牛頓法需要推導公式教材-應用MATLAB於數(shù)值計算u用兩次迭代解構造出的有限差分近似，替代牛頓法中的求導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)點u不需要顯式計算f(x) ，而有類似收斂的性質。en+1=O(enen-1)，是超線性收斂

9、，事實上en+1=O(en)u=(1+5)/2111()()(),nnnnnnnnnf xf xf xsxxxxs教材-應用MATLAB於數(shù)值計算u既然割線法使用前兩個近似解得到下一個解，那麼何既然割線法使用前兩個近似解得到下一個解，那麼何不利用前三個近似解呢不利用前三個近似解呢? ?u作法：作法：u用三個點用三個點a,b,ca,b,c及其函數(shù)值及其函數(shù)值f(a), f(b), f(c)f(a), f(b), f(c)，組一，組一個個2 2次曲線次曲線u對數(shù)據(jù)進行插值，然後另其與對數(shù)據(jù)進行插值，然後另其與x x軸的交點為下一個迭代軸的交點為下一個迭代點點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缺點缺點u速度不穩(wěn)定，需要好好馴服速度不穩(wěn)定，需要好好馴服( (人工多次反復測試人工多次反復測試) )教材-應用MATLAB於數(shù)值計算u結合二分法的可靠性與割線法及逆二次差值法結合二分法的可

11、靠性與割線法及逆二次差值法(IQI)(IQI)算法的收斂速度算法的收斂速度u演算法：演算法：u選取初始值選取初始值a a和和b b，使得，使得f(a)f(a)和和f(b)f(b)的正負號正好相反的正負號正好相反u使用一步割線法，得到使用一步割線法，得到a a和和b b之間的一個值之間的一個值u重復下面的步驟，直到重復下面的步驟，直到|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; 教材-應用MATLAB於數(shù)值計算 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; 教材-應用MATLAB於數(shù)值計算 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 教材-應用MATLAB於數(shù)值計算 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教材-應用MATLAB於數(shù)值計算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ū)間教材-應用MATLAB於數(shù)值計算二分點IQI點割線點教材-應用MATLAB於數(shù)值計算1.給定一個函數(shù)F(x)和值，求使得F()=;2.給定對未知函數(shù)F(x)採樣得到的一些數(shù)據(jù)點(xk,yk)和值，求使得F()=;解：1.

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

17、(-1)75步u後，區(qū)間步長大致減小為原始長度的eps倍。u經(jīng)過開始的一些步後，通常有足夠的歷史訊息，給出區(qū)間內(nèi)3個不同的點，以及對應的函數(shù)值。u如果這3個點產(chǎn)生的拋物線的極小(大)值落在區(qū)間內(nèi)，那麼下一點通常選擇這個極小點，而不是黃金分割點。2131011(35) / 20.382 教材-應用MATLAB於數(shù)值計算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極大值教材-應用MATLAB於數(shù)值計算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