数值计算方程式求根篇.ppt

1,數值計算 方程式求根篇,陳鴻智,2,有關於本教材,數值計算是件很好玩的事，可以用很精簡的程式完成 並且可以將它用圖形的方式展示出來。 Matlab可以讓我們作到。,3,教材內容,數值計算介紹 試圖告訴你數值計算是件很棒的事,4,為啥需要數值計算?,工程需求 資料探索分析 “透過觀察讓科學大幅進展”牛頓 “數學就是要求嚴密、美感和結構” - 高斯,牛頓,高斯,5,George Forsythe 於1960年代後期在 Stanford 大學開創了數值方法的課程。 Forsythe, Malcolm和 Moler合寫了一本教材。 那時的課程都是以Fortran語言編寫的程式庫為主。 為解決線性代數課程中矩陣運算問題， Cleve Moler博士於1980年開發出一套簡單易用的矩陣運算工具， Matlab於是誕生。,MATLAB語言的發展,Cleve Moler,6,歷史故事Matlab(MATrix LABoratory),由 Cleve Moler 研發將EISPACK 2,000條Fortran語法程式碼 80 個函數 (現已超過 8000個函數 ) 1979: 在史丹佛與 Jack Little 相遇，開始將源碼轉換成C 1984: Mathworks 成立,7,Matlab與Mathworks現況,Matlab, Simulink, Stateflow Matlab 2009b 各種領域的運用 航太,國防, 電腦,通信, 測繪 Mathworks 公司 3,500 以上的大學使用,500,000以上的用戶,100以上國家使用 2008年收益: 3億美元以上. 2008員工: 2,000人以上 價錢: 2,000$ (商務版價), 100$ (學生版價),8,誰需要?,研發單位 老師 實驗室軟體設備之一 Matlab 著重在演算法而不是程式語言探討 學生 批次處理檔案 不必再需要另外學其他的語言! 給你一個測試概念的好環境 把概念很快轉換成程式碼, 然後轉成 C或Java 等 易懂的視覺圖,9,應具備的課程知識,微積分 常微分方程 矩陣 計算機程式編寫經驗,10,二分法,是求實數自變量、實值函數f(x)=0的解之較慢但很考靠的算法 f(x)是個連續函數，且在a,b區間內會改變正負號，則可以找到一個x*，使得f(x*)=0 MATLAB k = 0; while abs(b-a) eps*abs(b) x = (a + b)/2; if sign(f(x) = sign(f(b) b = x; else a = x; end k = k + 1; end,11,牛頓法,公式xn+1= xn-f(xn)/f(xn) MATLAB k = 0; while abs(x - xprev) eps*abs(x) xprev = x; x = x - f(x)/fprime(x) k = k + 1; end 以計算為 例 while abs(x - xprev) eps*abs(x) xprev = x; x = 0.5*(x + M/x); end 缺點 當f(x)不連續、無有界的一階、二皆導數，或者初始值不接近準確解時就無法使用本方法,無限循環的牛頓法,需要推導公式,12,割線法,用兩次迭代解構造出的有限差分近似，替代牛頓法中的求導 公式 MATLAB while abs(b - a) eps*abs(b) c = a; a=b; b = b + (b-c)/f(c)/f(b)-1); k = k + 1; end 優點 不需要顯式計算f(x) ，而有類似收斂的性質。en+1=O(enen-1)，是超線性收斂，事實上en+1=O(en) =(1+5)/2,13,逆二次插值,既然割線法使用前兩個近似解得到下一個解，那麼何不利用前三個近似解呢? 作法： 用三個點a,b,c及其函數值f(a), f(b), f(c)，組一個2次曲線 對數據進行插值，然後另其與x軸的交點為下一個迭代點 MATLAB k=0; while abs(c - b) eps*abs(c) x = ployinterp(f(a),f(b),fc,a,b,c,0) a=b; b = c c=x; k = k + 1; end 缺點 速度不穩定，需要好好馴服(人工多次反復測試),14,Zeroin算法(一種綜合算法,結合它法的優點),結合二分法的可靠性與割線法及逆二次差值法(IQI)算法的收斂速度 演算法： 選取初始值a和b，使得f(a)和f(b)的正負號正好相反 使用一步割線法，得到a和b之間的一個值 重復下面的步驟，直到|b -a|b|或者f(b) =0 重新排列a、b和c(可能要經過兩輪)，使得 f(a)和f(b)的正負號相反 |f(b) |f(a) | c的值為上一個b值 如果c a，執行IQI算法中的一步疊代 如果c= a ，執行割線法中的一步 如果執行IQI算法或割線法得到的近似解在區間a,b內，接受這個解為c 如果這個近似解不在區間a,b內，執行一步二分法得到c 。,15,Zeroin算法,function b = fzerotx(F,ab,varargin) if ischar(F) ,16,Zeroin算法, while fb = 0 if sign(fa) = sign(fb) a = c; fa = fc; d = b - c; e = d; end if abs(fa) c = b; b = 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) d = m; e = m; else s = fb/fc;,17,Zeroin算法, 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 - 1.0)*(r - 1.0)*(s - 1.0); end; if p 0, q = -q; else p = -p; end; if (2.0*p e = d; d = p/q; else d = m; e = m; end; end,18,Zeroin算法, 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,19,Zeroin算法,MATLAB for n = 1:10 z(n) = fzerotx(besselj(0,x),(n-1) n*pi); end x = 0:pi/50:10*pi; y = besselj(0,x); plot(z,zeros(1,10),o,x,y,-) line(0 10*pi,0 0,color,black) axis(0 10*pi -0.5 1.0),輸入:基本函數、搜索區間,20,fzerogui-看圖識zeroin演算法,二分點,IQI點,割線點,21,尋找函數為某個值的解和反向插值,1.給定一個函數F(x)和值，求使得F()=; 2.給定對未知函數F(x)採樣得到的一些數據點(xk,yk)和值，求使得F()=; 解： 1. f(x)=F(x)-=0，可以求得x= 2. f(x)=P(x)-， P(x)是某個插值函數(可以用pchiptx(xk,yk,x)或splinetc(xk,yk,x)得到P(x) ，其缺點為：計算代價大，需要反復計算插值。 使用反復插值法，用的是pchip和spline算法，將xk和yk角色互調 yk需具單調性，或者至少yk的某個包含目標值的子集單調。 生出另一個分段多項式Q(yk)=xk。y=，=Q(y)。,22,最優化和fmintx,f(x),和它所在的區間a,b ，求x使得f(x)為極小(或極大) 作法: 黃金分割搜索+拋物線法 區間步長會以-10.618的比例 縮小，經過-52/log2(-1)75步 後，區間步長大致減小為原始長度的eps倍。 經過開始的一些步後，通常有足夠的歷史訊息，給出區間內3個不同的點，以及對應的函數值。 如果這3個點產生的拋物線的極小(大)值落在區間內，那麼下一點通常選擇這個極小點，而不是黃金分割點。,23,最優化和fmintx,求 的極大值 F = inline(-humps(x); fmintx(F,-1,2,1.e-4) Step x f(x) init: 0.1458980337 -25.2748253202 gold: 0.8541019662 -20.9035150009 gold: -0.2917960675 2.5391843579 para: 0.4492755129 -29.0885282699 para: 0.4333426114 -33.8762343193 para: 0.3033578448 -96.4127439649 gold: 0.2432135488 -71.7375588319 para: 0.3170404333 -93.8108500149 para: 0.2985083078 -96.4666018623 para: 0.3003583547 -96.5014055840 para: 0.3003763623 -96.5014085540 para: 0.3003756221 -96.5014085603,極大值,24,進一步閱讀,R. P. Brent, Algorithms for Minimization Without Derivatives, Prentice Hall,Englewood Cliffs, 1973. G. Dahlquest and A. Bjorck, Numerical Methods, Prentice Hall, Englewood Cliffs, 1974. T. J. Dekker, Finding a Zero by Means of Successive Linear Interpolation, Constructive Aspects of the