1、會計學1Exp姜啟源數學建模珍藏姜啟源數學建模珍藏計算機會“算”嗎？靠得住嗎？例：把例：把4開開n次方，再平方次方，再平方n次，結果是次，結果是4？存在誤差？存在誤差？英國著名數值分析學家英國著名數值分析學家 Higham (1998): Higham (1998): Can you count on computers?Can you count on computers?精確計算：精確計算：解析結果解析結果 (Analytical)近似近似計計算：算：數值數值結果結果(Numerical)?422n=55左右：結果變成左右：結果變成1計計算功效算功效=計算工具計算工具*計算方法計算方法(算

2、法算法)浮點運算：舍入誤差浮點運算：舍入誤差第1頁/共42頁實驗3的基本內容3.3.數值積分的數值積分的梯形公式、辛普森公式和梯形公式、辛普森公式和高斯公式。高斯公式。1.1.插值的基本原理；插值的基本原理； 三種插值方法：拉格朗日插三種插值方法：拉格朗日插 值，分段線值，分段線性性 插值，三次樣條插值。插值，三次樣條插值。2.2.插值的插值的 MATLAB 實現實現及插值的應用及插值的應用。4.4.數值積分的數值積分的 MATLAB 實現實現及數值積分的應用及數值積分的應用。第2頁/共42頁什么是插值什么是插值(Interpolation)？從查函數表說？從查函數表說起起查查 函函 數數 表

3、表xtdtex2221)(x0121.0 0.8413 0.8438 0.84611.1 0.8643 0.8665 0.86861.2 0.8849 0.8869 0.8888標準正態分布函數表標準正態分布函數表求求 (1.114) (1.114)=0.8665 (0.8686 0.8665) 0.4=0.8673插值插值插值在圖像處理插值在圖像處理/數控加工數控加工/外觀設計等領域有重要應用外觀設計等領域有重要應用第3頁/共42頁插值的基本原理插值的基本原理插值問題的提法插值問題的提法已知已知 n+1n+1個節點個節點, 1 , 0(),(njyxjj其中其中jx互不相同，不妨設互不相同，

4、不妨設),10bxxxan求任一插值點求任一插值點)(*jxx 處的插值處的插值.*y0 x1xnx0y1y節點可視為由節點可視為由)(xgy 產生產生,g表達式復雜表達式復雜,甚至無表達甚至無表達式式*x*y第4頁/共42頁0 x1xnx0y1y求解插值問題的基本思路求解插值問題的基本思路構造一個構造一個( (相對簡單的相對簡單的) )函數函數),(xfy 通過全部節點通過全部節點, ,即即), 1 ,0()(njyxfjj再用再用)(xf計算插值，計算插值，即即).(*xfy *x*y插值的插值的基本原理基本原理第5頁/共42頁1.1.拉格朗日拉格朗日(Lagrange)(Lagrange

5、)多項式插值多項式插值1.0 1.0 插值多項式插值多項式) 1 ()(0111axaxaxaxLnnnnnnnnnnnnnyyYaaAxxxxX001100,11在什么條件下）(0)det(X), 1 , 0()(njyxLjjn)2(YXA 求ia三種插值三種插值方法方法有唯一解)2(第6頁/共42頁1.1 1.1 拉格朗日插值多項拉格朗日插值多項式式nixxxxxxxxxxxxxxxxxlniiiiiiniii1 , 0,)()()()()()()(110110)3()()(0 xlyxLiniinjjnjiyxLjijixl)(,0, 1)(則若又又（2）有唯一解，故有唯一解，故（3）

6、與與（1）相同相同。 基函數基函數( )ilx) 1 ()(0111axaxaxaxLnnnnn) 2(YXA三種插值三種插值方法方法第7頁/共42頁),(),()!1()()()()(0)1(baxxngxLxgxRnjjnnn1)1()(nnMg減小（粗略地看）如何使誤差)(xRn平緩gjxx 接近njjnnxxnMxR01)!1()(三種插值三種插值方法方法1.2 1.2 誤差估計誤差估計增加n第8頁/共42頁1.3 1.3 拉格朗日插值多項式的振蕩拉格朗日插值多項式的振蕩?)(?)(xRxLnnn55,11)(2xxxg63. 363. 3),()(limxxgxLnnRunge現現象

7、象取n=2,4,6,8,10,計算Ln(x), 畫出圖形-505-1.5-1-0.500.511.52y=1/(1+x2)n=2n=4n=6n=8n=10三種插值三種插值方法方法Runge.m第9頁/共42頁2.2.分段線性插值分段線性插值xjxj-1xj+1x0 xn其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxI計算量與計算量與n n無關無關; ;n n越大，誤差越小越大，誤差越小. .nnnxxxxgxI0),()(lim三種插值方法三種插值方法第10頁/共42頁機翼下輪廓線3. 3. 三次樣條插值三次樣條插值樣條函數的由來樣

8、條函數的由來飛機、船體、汽車外形等的放樣（設計飛機、船體、汽車外形等的放樣（設計）細木條：樣條細木條：樣條第11頁/共42頁3. 3. 三次樣條插值三次樣條插值, 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs數學樣條（數學樣條（spline）iiiidcban,4 個待定系數3）) 1, 1()()()()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii3）2），3）共 4n-2個方程三種插值方法三種插值方法第12頁/共42頁自然邊界條件）(0)()()40

9、 nxSxS)(,)4)3)2xSdcbaiiii三次樣條插值確定三次樣條插值確定4 4n n個系數需增加個系數需增加 2 2個條件個條件思考1）自然邊界條件的幾何意義是什么？2）樣條插值為什么普遍用3次多項式，而不是2或4次？三次樣條插值三次樣條插值).()(limxgxSn第13頁/共42頁三種插值方法小結三種插值方法小結 拉格朗日插值（高次多項式插值）：拉格朗日插值（高次多項式插值）：曲線光滑；誤差估計有表達式；收斂性不能保證。曲線光滑；誤差估計有表達式；收斂性不能保證。用于理論分析，實際意義不大用于理論分析，實際意義不大。 分段線性和三次樣條插值（低次多項式插值）：分段線性和三次樣條插

10、值（低次多項式插值）：曲線不光滑（三次樣條插值已大有改進）；誤差估曲線不光滑（三次樣條插值已大有改進）；誤差估計較難（對三次樣條插值）；收斂性有保證。計較難（對三次樣條插值）；收斂性有保證。簡單實用，應用廣泛簡單實用，應用廣泛。 其他：其他：Hermite插值、分段三次插值、二維插值等插值、分段三次插值、二維插值等根據需要，根據需要，各取所需各取所需。第14頁/共42頁1. 1. 拉格朗日插值拉格朗日插值: :自編程序自編程序, ,如名為如名為 lagr.m 的的M文件，文件， 第一行為第一行為 function y=lagr(x0,y0,x) 輸入輸入: :節點節點x0,y0, 插值點插值點

11、x ( (均為均為數組，長度自定義數組，長度自定義) )）；）； 輸出輸出: :插值插值y ( (與與x同長度數組同長度數組) )）。）。 應用時輸入應用時輸入x0,y0,x后后, ,運行運行 y=lagr(x0,y0,x)2. 2. 分段線性插值分段線性插值: :已有程序已有程序 y=interp1(x0,y0,x) y=interp1(x0,y0,x,linear)3. 3. 三次樣條插值三次樣條插值: :已有程序已有程序 y=interp1(x0,y0,x,spline) 或或 y=spline(x0,y0,x)用MATLAB作插值計算注：注：MATLAB有樣條工具箱（有樣條工具箱（Sp

12、line Toolbox）第15頁/共42頁用MATLAB作插值計算55,11)(2xxxg為例，作三種插值的比較為例，作三種插值的比較以以 0 1.0000 1.0000 1.0000 1.0000 0.5000 0.8000 0.8434 0.7500 0.8205 1.0000 0.5000 0.5000 0.5000 0.5000 1.5000 0.3077 0.2353 0.3500 0.2973 2.0000 0.2000 0.2000 0.2000 0.2000 2.5000 0.1379 0.2538 0.1500 0.1401 3.0000 0.1000 0.1000 0.1

13、000 0.1000 3.5000 0.0755 -0.2262 0.0794 0.0745 4.0000 0.0588 0.0588 0.0588 0.0588 4.5000 0.0471 1.5787 0.0486 0.0484 5.0000 0.0385 0.0385 0.0385 0.0385 x y y1 y2 y3 用用n=11個節點，個節點，m=21個插值點，三種方法作插值，畫圖。個插值點，三種方法作插值，畫圖。chazhi1第16頁/共42頁插值的應用加工時需要加工時需要x每改變每改變0.05時的時的y值值chazhi2圖1 零件的輪廓線 （x間隔0.2）表1 x間隔0.2的加

14、工坐標x,y（圖1右半部的數據）數控機床加工零件數控機床加工零件 0.0,5.00 0.2,4.710.4,4.31 0.6,3.68 0.8,3.051.0,2.50 1.2,2.051.4,1.69 1.6,1.40 1.8,1.182.0,1.00 2.2,0.862.4,0.74 2.6,0.64 模型模型 將圖1逆時針方向轉90度，輪廓線上下對稱，只需對上半部計算一個函數在插值點的值。 圖2 逆時針方向轉90度的結果-5-4-3-2-101234500.511.522.533.544.55uv令v=x, u= -y 第17頁/共42頁為什么要作數值積分為什么要作數值積分 許多函數許多

15、函數“積不出來積不出來”, ,只能用數值方法，如只能用數值方法，如dxxxdxebabaxsin,22 積分是重要的數學工具，是微分方程、概率積分是重要的數學工具，是微分方程、概率論等的基礎；在實際問題中有直接應用。論等的基礎；在實際問題中有直接應用。 對于用離散數據或者圖形表示的函數對于用離散數據或者圖形表示的函數，(可以先做插值然后積分；或者直接利用點做數值積分)計算積分只有求助于數值方法。計算積分只有求助于數值方法。數值積分數值積分第18頁/共42頁nabfIIdxxfInkknnnba)(,lim)(1數數 值值 積積 分分 的的 基基 本本 思思 路路回回 憶憶 定定 積積 分分 的

16、的 定定 義義各種數值積分方法研究的是各種數值積分方法研究的是k),(ba如何取值，區間如何取值，區間如何劃分，如何劃分，使得既能保證一定精度，計算量又小。使得既能保證一定精度，計算量又小。n n充分大時充分大時I In n就是就是I I的數值積分的數值積分（計算功效：算得準，算得快）（計算功效：算得準，算得快）第19頁/共42頁1.1.從矩形公式到梯形公式數值積分數值積分yy=f(x)xbao)1(10nkknfhL)(,10kknkxffnabhbxxxxa)2(1nkknfhRnnRL,平均，得到梯形公式) 3()(2011nnkknffhfhTxk+1xkxk-1fk第20頁/共42頁

17、2.2.辛普森辛普森(Simpson)(Simpson)公式公式（拋物線公式）（拋物線公式） 梯形公式相當于用分段線性插值函數分段線性插值函數代替)(xf每段要用相鄰兩小區間兩小區間端點的三個函數值端點的三個函數值拋物線拋物線公式公式提高精度提高精度分段二次插值函數分段二次插值函數2221212222(,),(,),(,)0,1, ,1kkkkkkxfxfxfkm數值積分數值積分yy=f(x)xbaox2kf2kx2k+1x2k+2f2k+1f2k+2區間數必須為偶數區間數必須為偶數mn2第21頁/共42頁) 4(2),24(3112101220mabhffffhSmkkmkkmm 對對k求和

18、求和(共共m段段) )，得（復合），得（復合）辛普森公式辛普森公式：)4(3)(22122222kkkxxkfffhdxxskk二次插值函數sk(x)構造用),(),(),(2222121222kkkkkkfxfxfx2.2.辛普森辛普森(Simpson)(Simpson)公式（拋物線公式）公式（拋物線公式）第22頁/共42頁bannnTdxxfTITfR)(),(梯形公式在每小段上是用梯形公式在每小段上是用線性插值函數線性插值函數T T( (x) )代替代替 f( (x) )11()( )( )()(),(,)2kkkkkkff xT xxxxxxx x （拉格朗日插值余項）梯形公式梯形公式

19、的誤差估計的誤差估計)(2011nnkknffhfhTbadxxf)()(12)(2)()()(3111kxxkkkxxfhdxxxxxfdxxTxfkkkk 因為：因為：(x-xk)(x-xk+1)在在(xk,xk+1)不變號，所以：不變號，所以：第23頁/共42頁)5()(12|),(|22abMhTfRn梯形公式梯形公式Tn的的誤差是誤差是h2階的階的),(, )(max2baxxfM 估計估計habn因為 103)(12|),(|nkknfhTfR 梯形公梯形公式的誤差式的誤差)( )( (121)(121 2afbfdxxfhTIban 103)(12)(nkkbannfhTdxxf

20、TI) 5()()(122afbfhTIn第24頁/共42頁同理可得：同理可得：) 6()(180| ),(|44abMhSfRn其中其中),(,)(max)4(4baxxfM辛普森公式辛普森公式Sn的誤差是的誤差是h4階的階的。辛普森公式的誤差估辛普森公式的誤差估計計第25頁/共42頁梯形公式和辛普森公式的收斂性若對若對I某個數值積分某個數值積分In有有chIIpnnlim（非零常數）（非零常數）則稱則稱 In是是 p 階收斂的階收斂的。梯形公式梯形公式 2 2 階收斂，辛普森公式階收斂，辛普森公式 4 4 階收斂。階收斂。c=0: 至少至少p階收斂（超階收斂（超p階收斂）階收斂）第26頁/

21、共42頁積分步長的自動選取積分步長的自動選取選定數值積分公式后，如何確定步長選定數值積分公式后，如何確定步長h以滿足給定的誤差以滿足給定的誤差 )()(122afbfhTIn梯形公式)(412nnTITInnTT2用二分法只要用二分法只要其中其中fk+1/2是原是原分點分點xk,xk+1的中點的中點（記記xk+1/2）的函數值的函數值1021222nkknnfhTT且且T T2n2n可在可在T Tn n基礎上計算基礎上計算)(3122nnnTTTInTI2)2(2nnhh第27頁/共42頁高斯高斯(Gauss)(Gauss)求積公式求積公式矩形公式矩形公式(1)、(2)梯形公式梯形公式(3)辛

22、普森公式辛普森公式(4)A Ak k是與是與f f無關的常數無關的常數代數代數精度精度設設,)(kxxf用用(7)計算計算,)(badxxfI若對于若對于mk,1 ,0都有都有, IIn而當而當, 1IImkn則則稱稱In的代數精度為的代數精度為m.)7()(1nkkknxfAINewton-Cotes方法方法第28頁/共42頁梯形公式的代數精度（考察梯形公式的代數精度（考察T1）k=1f(x)=x222abxdxIba2)()()(21baabbfafhT3332abdxxIba2)(221baabTk=2f(x)=x2IT1IT 1梯形公式的代數精度為梯形公式的代數精度為1辛普森公式的代數

23、精度為辛普森公式的代數精度為3第29頁/共42頁高斯公式的思路高斯公式的思路取消對節點的限制，按照代數精度最大取消對節點的限制，按照代數精度最大的原則，同時確定節點的原則，同時確定節點xk和系數和系數Ak構造求積公式構造求積公式)()(22112xfAxfAG對于對于11)(dxxfI使使G G2 2的代數精度為的代數精度為3 332, 1)(xxxxf)()()(221111xfAxfAdxxf確定確定2121,AAxx第30頁/共42頁03/202322311222211221121xAxAxAxAxAxAAA將將f(x)f(x)代入計算得代入計算得1,3/1,3/12121AAxx)3/

24、1 ()3/1(2ffG用用n個節點，個節點，Gn的代數精度可達的代數精度可達2n-1, 但是需解但是需解復雜的非線性方程組，實用價值不大。復雜的非線性方程組，實用價值不大。第31頁/共42頁常 用 的 高 斯 公 式將將( (a,ba,b) )分小，把小區間變換為分小，把小區間變換為(-1(-1，1), 1), 再用再用G G2 2mkkkbazfzfhdxxf1)2()1 ()()(2)(1)(2)11,222 32 3kkkkkkxxxxhhzz伸縮變換公式mkkhaxmabhk, 1, 0,/ )(代數精度為代數精度為3節點加密時，原計算信息無法利用節點加密時，原計算信息無法利用第32

25、頁/共42頁思路思路：將積分區間分小，在小區間上用：將積分區間分小，在小區間上用n不太不太 大大的的 。而在節點加密一倍時能夠利用原節點的。而在節點加密一倍時能夠利用原節點的函數值，可以把區間的端點作為固定節點。函數值，可以把區間的端點作為固定節點。改進的高斯公式nG)()()(121bfAxfAafAGnnkkknGauss-Lobatto求積公式求積公式 其中a, b為小區間的端點，nnAAxx,112為2n-2個參數，代數精度可達到代數精度可達到2n-3注意：實際計算中一般采用自適應方法確定步長注意：實際計算中一般采用自適應方法確定步長第33頁/共42頁用用MATLAB 作數值積分作數值

26、積分10nkknfhLnkknfhR1矩形矩形公式公式Sum(x)輸入數組x(即fk),輸出x的和(數)cumsum(x)輸入數組x,輸出x的依次累加和(數組)梯形梯形公式公式)(2011nnkknffhfhTTrapezoidal (梯形)trapz(x)輸入數組x,輸出按梯形公式x的積分(單位步長)trapz(x,y)輸入同長度數組 x,y,輸出按梯形公式y對x的積分(步長不一定相等)第34頁/共42頁用用MATLAB 作數值積分作數值積分mabhffffhSmkkmkkmn2),24(3112101220辛普森公式辛普森公式(quadrature: 求積) quad(fun,a,b,to

27、l,trace)I,fn=quad()用自適應辛普森公式計算tol為絕對誤差，缺省時為10-6Gauss-Lobatto公公式式)()()(121bfAxfAafAGnnkkknquadl(fun,a,b,tol,trace)I,fn=quadl()用自適應Gauss-Lobatto公式計算 tol為絕對誤差，缺省時為10-6注意：注意：fun.m中應以自變量為矩陣的形式輸入中應以自變量為矩陣的形式輸入(點運算點運算)第35頁/共42頁矩形域上計算二重積分的命令：矩形域上計算二重積分的命令：dblquad(fun,xmin,xmax,ymin,ymax,tol)廣義積分、二重和三重積分長方體上計算三重積分的命令：長方體上計算三重積分的命令：triplequad(fun,xmin,xmax,ymin,ymax, zmin,zmax,tol)注：注：fun是被積函數，本身可以有自己的參數是被積函數，本身可以有自己的參數廣義積分：廣義積分：通過分析和控制誤差，轉換成普通積分通過分析和控制誤差，轉換成普通積分quadv(fun,a,b,tol,trace)向量值積分：向量值積分：第36頁/共42頁用用MATLAB 作數值積分作數值積分例例. 計算計算4011s i nd xx1 1）矩形公式和梯形公