1、大學數學實驗大學數學實驗Mathematical Experiments 實驗實驗3 3 插值與數值積分插值與數值積分計算機會“算”嗎？靠得住嗎？例：把例：把4開開n次方，再平方次方，再平方n次，結果是次，結果是4？存在誤差？存在誤差？英國著名數值分析學家英國著名數值分析學家 Higham (1998): Higham (1998): Can you count on computers?Can you count on computers?精確計算：精確計算：解析結果解析結果 (Analytical)近似計算：近似計算：數值結果數值結果(Numerical)?422n=55左右：結果變成左右

2、：結果變成1計算功效計算功效=計算工具計算工具*計算方法計算方法(算法算法)浮點運算：舍入誤差浮點運算：舍入誤差實驗3的基本內容3.3.數值積分的梯形公式、辛普森公式和高斯公式。數值積分的梯形公式、辛普森公式和高斯公式。1.1.插值的基本原理；插值的基本原理； 三種插值方法：拉格朗日插三種插值方法：拉格朗日插 值，分段線性值，分段線性 插值，三次樣條插值。插值，三次樣條插值。2.2.插值的插值的 MATLAB 實現及插值的應用。實現及插值的應用。4.4.數值積分的數值積分的 MATLAB 實現及數值積分的應用。實現及數值積分的應用。什么是插值什么是插值(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插值插值插值在圖像處理插值在圖像處理/數控加工數控加工/外觀設計等領域有重要應用外觀設計等領域有重要應用插值的基本原理插值的基本原理插值問題的提法插值問題的提法已知已知 n+1n+1個節點個節點, 1 , 0(),(njyxjj其中其中jx互

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

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

6、基函數基函數( )ilx) 1 ()(0111axaxaxaxLnnnnn)2(YXA三種插值三種插值方法方法9思考1.對于n+1個節點，若用次數大于或小于n的多項式作插值，結果如何？2.由n+1個節點得到的Ln(x)的次數會不會小于n?試就n=2的情況加以說明。3.若g(x)為m次多項式， ,問mn( )( )nL xg x與關系如何？10 對于給定的n+1個點 ，對應于它們的次數不超過n的拉格朗日多項式Ln(x)只有一個。如果計入次數更高的多項式，則有無窮個，因為所有與Ln(x)相差 的多項式都滿足條件。01()().()nxxxxxx11練習( )sin(/2)0,1n+1n=1,2,3

7、 lagrange yg xx用函數 在區間 上等距地產生 個節點，當 時計算 插值多項式。),(),()!1()()()()(0)1(baxxngxLxgxRnjjnnn1)1()(nnMg減小（粗略地看）如何使誤差)(xRn平緩gjxx 接近njjnnxxnMxR01)!1()(三種插值三種插值方法方法1.2 1.2 誤差估計誤差估計增加n1.3 1.3 拉格朗日插值多項式的振蕩拉格朗日插值多項式的振蕩?)(?)(xRxLnnn55,11)(2xxxg63. 363. 3),()(limxxgxLnnRunge現象現象取n=2,4,6,8,10,計算Ln(x), 畫出圖形-505-1.5-

8、1-0.500.511.52y=1/(1+x2)n=2n=4n=6n=8n=10三種插值三種插值方法方法Runge.m2.2.分段線性插值分段線性插值xjxj-1xj+1x0 xn其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxI計算量與計算量與n n無關無關; ;n n越大，誤差越小越大，誤差越小. .nnnxxxxgxI0),()(lim三種插值三種插值方法方法MATLABxch11，xch12，xch13，xch1466,11)(2xxxg例用分段線性插值法求插值,并觀察插值誤差.1.在-6,6中平均選取5個點作插值(xch11

9、)4.在-6,6中平均選取41個點作插值(xch14)2.在-6,6中平均選取11個點作插值(xch12)3.在-6,6中平均選取21個點作插值(xch13)比分段線性插值更光滑xyxi-1 xiab 在數學上，光滑程度的定量描述是：函數(曲線)的k階導數存在且連續，則稱該曲線具有k階光滑性 光滑性的階次越高，則越光滑。是否存在較低次的分段多項式達到較高階光滑性的方法？三次樣條插值就是一個很好的例子。三次樣條插值機翼下輪廓線3. 3. 三次樣條插值三次樣條插值樣條函數的由來樣條函數的由來飛機、船體、汽車外形等的放樣（設計）飛機、船體、汽車外形等的放樣（設計）細木條：樣條細木條：樣條3. 3.

10、三次樣條插值三次樣條插值, 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs數學樣條（數學樣條（spline）iiiidcban,4 個待定系數3）111( )( ),( )( )( )( )(1,21)iiiiiiiiiiiis xsxs xsxs xsxin，.2），3）共 4n-2個方程三種插值三種插值方法方法3）自然邊界條件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiii三次樣條插值確定三次樣條插值確定4 4n n個系數需增加個系數需增加 2 2個條件個條

11、件三次樣條三次樣條插值插值).()(limxgxSn思考1）自然邊界條件的幾何意義是什么？自然樣條是通過所有數據點的插值函數中，曲率最小的唯一函數。因此，自然三次樣條是插值所有數據點的最光滑的函數。2）樣條插值為什么普遍用3次多項式，而不是2或4次？樣條函數不一定必須是逐段三次多項式，也可以是逐段的簡單函數，且連接點保持足夠光滑。但因三次多項式計算簡單，且能滿足一般實際問題的需要，故三次樣條函數用的最多。例66,11)(2xxxg用三次樣條插值選取11個基點計算插值(ych)To MATLAB ych(larg1)三種插值方法小結三種插值方法小結 拉格朗日插值（高次多項式插值）：拉格朗日插值（

12、高次多項式插值）：曲線光滑；誤差估計有表達式；收斂性不能保證。曲線光滑；誤差估計有表達式；收斂性不能保證。用于理論分析，實際意義不大用于理論分析，實際意義不大。 分段線性和三次樣條插值（低次多項式插值）：分段線性和三次樣條插值（低次多項式插值）：曲線不光滑（三次樣條插值已大有改進）；誤差估曲線不光滑（三次樣條插值已大有改進）；誤差估計較難（對三次樣條插值）；收斂性有保證。計較難（對三次樣條插值）；收斂性有保證。簡單實用，應用廣泛簡單實用，應用廣泛。 其他：其他：Hermite插值、分段三次插值、二維插值等插值、分段三次插值、二維插值等根據需要，各取所需根據需要，各取所需。1. 1. 拉格朗日插

13、值拉格朗日插值: :自編程序自編程序, ,如名為如名為 lagr.m 的的M文件，文件， 第一行為第一行為 function y=lagr(x0,y0,x) 輸入輸入: :節點節點x0,y0, 插值點插值點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. 三次樣條插

14、值三次樣條插值: :已有程序已有程序 y=interp1(x0,y0,x,spline) 或或 y=spline(x0,y0,x)用MATLAB作插值計算注：注：MATLAB有樣條工具箱（有樣條工具箱（Spline Toolbox）用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

15、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.1000 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個插值點，個插值點，三種方法作三種方法作插值，

16、畫圖。插值，畫圖。chazhi1插值的應用加工時需要加工時需要x每每改變改變0.05時的時的y值值chazhi2圖1 零件的輪廓線 （x間隔0.2）表1 x間隔0.2的加工坐標x,y（圖1右半部的數據）數控機床加工零件數控機床加工零件 模型模型 將圖1逆時針方向轉90度，輪廓線上下對稱，只需對上半部計算一個函數在插值點的值。 圖2 逆時針方向轉90度的結果-5-4-3-2-101234500.511.522.533.544.55uv令v=x, u= -y 例：從1點到12點的11小時內，每隔1小時測量一次溫度，測得的溫度的數值依次為：5，8，9，15，25，29，31，30，22，25，27，

17、24試估計每隔1/10小時的溫度值。To MATLAB (temp)hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline); (直接輸出數據將是很多的)plot(hours,temps,+,h,t,hours,temps,r:) %作圖xlabel(Hour),ylabel(Degrees Celsius)xy機翼下輪廓線例 已知飛機下輪廓線上數據如下，求x每改變0.1時的y值。To MATLAB(plane)二維插值的定義xyO第一種（網格節點）： 已知 mn個節

18、點 (,) (1,2,.,;1,2, )ijijxyzim jn其中jiyx ,互不相同，不妨設bxxxam 21dyyycn 21 構造一個二元函數),(yxfz 通過全部已知節點,即再用),(yxf計算插值，即).,(*yxfz (,)(0,1,;0,1,)ijijfxyzimjn第二種（散亂節點）：yxO已知n個節點),.,2 , 1(),(nizyxiii 其中),(iiyx互不相同， 構造一個二元函數),(yxfz 通過全部已知節點,即),1 ,0(),(nizyxfiii 再用),(yxf計算插值，即).,(*yxfz 注意：最鄰近插值一般不連續具有連續性的最簡單的插值是分片線性插

19、值最鄰近插值xy(x1, y1)(x1, y2)(x2, y1)(x2, y2)O 二維或高維情形的最鄰近插值，與被插值點最鄰近的節點的函數值即為所求 將四個插值點（矩形的四個頂點）處的函數值依次簡記為： 分片線性插值xy(xi, yj)(xi, yj+1)(xi+1, yj)(xi+1, yj+1)Of (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f4插值函數為：11()jjijiiyyyxxyxx12132( , )()()()()ijf x yfffxxffyy 第二片(上三角形區域)：(x, y)滿足11()jj

20、iiiiyyyxxyxx插值函數為：14134( , )()()()()jif x yfffyyffxx注意：(x, y)當然應該是在插值節點所形成的矩形區域內顯然，分片線性插值函數是連續的；分兩片的函數表達式如下：第一片(下三角形區域)： (x, y)滿足 雙線性插值是一片一片的空間二次曲面構成雙線性插值函數的形式如下：( , )()()f x yaxb cyd其中有四個待定系數，利用該函數在矩形的四個頂點（插值節點）的函數值，得到四個代數方程，正好確定四個系數雙線性插值xy(x1, y1)(x1, y2)(x2, y1)(x2, y2)O 要求x0,y0單調；x，y可取為矩陣，或x取行向量

21、，y取為列向量，x,y的值分別不能超出x0,y0的范圍z=interp2(x0,y0,z0,x,y,method)被插值點插值方法用MATLAB作網格節點數據的插值插值節點被插值點的函數值nearest 最鄰近插值；linear 雙線性插值；cubic 雙三次插值；缺省時 雙線性插值.例：測得平板表面35網格點處的溫度分別為： 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 試作出平板表面的溫度分布曲面z=f(x,y)的圖形輸入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 8

22、5 86;mesh(x,y,temps)1.先在三維坐標畫出原始數據，畫出粗糙的溫度分布曲線圖.2以平滑數據,在 x、y方向上每隔0.2個單位的地方進行插值.再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖. To MATLAB (wendu) 通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進行比較To MATLAB (moutain)返回 插值函數griddata格式為: cz =griddata（x，y，z，cx，cy，method）用MA

23、TLAB作散點數據的插值計算 要求cx取行向量，cy取為列向量被插值點插值方法插值節點被插值點的函數值nearest最鄰近插值linear 雙線性插值cubic 雙三次插值v4- MATLAB提供的插值方法缺省時, 雙線性插值 例 在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區域（75，200）（-50，150）里的哪些地方船要避免進入 2.75,20050,150. (1)hd 在矩形區域作二維插值三次插值法 .1 輸入插值基點數據To MATLAB hd14.作出水深小于5的海域范圍,即z=5的等高線. 3 作海底曲面圖432011高教社杯全國大學生數學建

24、模競賽A題A題題 城市表層土壤重金屬污染分析城市表層土壤重金屬污染分析將所考察的城區劃分為間距1公里左右的網格子區域，按照每平方公里1個采樣點對表層土（010 厘米深度）進行取樣、編號，并用GPS記錄采樣點的位置。應用專門儀器測試分析，獲得了每個樣本所含的多種化學元素的濃度數據。附件1列出了采樣點的位置、海拔高度及其所屬功能區等信息，附件2列出了8種主要重金屬元素在采樣點處的濃度。44現要求你們通過數學建模來完成以下任務：(1) 給出8種主要重金屬元素在該城區的空間分布， 并分析該城區內不同區域重金屬的污染程度。(2) 通過數據分析，說明重金屬污染的主要原因。(3) 分析重金屬污染物的傳播特征

25、，由此建立模型， 確定污染源的位置。(4) 分析你所建立模型的優缺點，為更好地研究城市 地質環境的演變模式，還應收集什么信息？有了 這些信息，如何建立模型解決問題？nabfIIdxxfInkknnnba)(,lim)(1數數 值值 積積 分分 的的 基基 本本 思思 路路回回 憶憶 定定 積積 分分 的的 定定 義義各種數值積分方法研究的是各種數值積分方法研究的是k),(ba如何取值，區間如何取值，區間如何劃分，如何劃分，使得既能保證一定精度，計算量又小。使得既能保證一定精度，計算量又小。n n充分大時充分大時I In n就是就是I I的數值積分的數值積分（計算功效：算得準，算得快）（計算功效

26、：算得準，算得快）為什么要作數值積分為什么要作數值積分 許多函數許多函數“積不出來積不出來”, ,只能用數值方法，如只能用數值方法，如dxxxdxebabaxsin,22 積分是重要的數學工具，是微分方程、概率積分是重要的數學工具，是微分方程、概率論等的基礎；在實際問題中有直接應用。論等的基礎；在實際問題中有直接應用。 對于用離散數據或者圖形表示的函數對于用離散數據或者圖形表示的函數，計算積分只有求助于數值方法。計算積分只有求助于數值方法。數值數值積分積分1.1.從矩形公式到梯形公式數值積分數值積分yy=f(x)xbao)1(10nkknfhL)(,10kknkxffnabhbxxxxa)2(

27、1nkknfhRnnRL,平均，得到梯形公式) 3()(2011nnkknffhfhTxk+1xkxk-1fk2.2.辛普森辛普森(Simpson)(Simpson)公式公式（拋物線公式）（拋物線公式） 梯形公式相當于用分段線性插值函數分段線性插值函數代替)(xf每段要用相鄰兩小區間兩小區間端點的三個函數值端點的三個函數值拋物線拋物線公式公式提高精度提高精度分段二次插值函數分段二次插值函數2221212222(,),(,),(,)0,1, ,1kkkkkkxfxfxfkm數值積分數值積分yy=f(x)xbaox2kf2kx2k+1x2k+2f2k+1f2k+2區間數必須為偶數區間數必須為偶數m

28、n2) 4(2),24(3112101220mabhffffhSmkkmkkmm 對對k求和求和(共共m段段) )，得（復合），得（復合）辛普森公式辛普森公式：)4(3)(22122222kkkxxkfffhdxxskk二次插值函數sk(x)構造用),(),(),(2222121222kkkkkkfxfxfx2.2.辛普森辛普森(Simpson)(Simpson)公式（拋物線公式）公式（拋物線公式）bannnTdxxfTITfR)(),(梯形公式在每小段上是用梯形公式在每小段上是用線性插值函數線性插值函數T T( (x) )代替代替 f( (x) ),(,),)(2)()()(11 kkkkk

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

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

31、科茲公式計算積分(精確值為 ).4 , 2 , 1n43096441. 06244267767. 0)15 . 0(25 . 0115 . 0dxx43093403. 0)175. 045 . 0(65 . 0115 . 0dxx43096407. 0) 17875. 03275. 012625. 0325 . 07 (905 . 0115 . 0dxx利用梯形公式利用拋物線公式利用樣條公式積分步長的自動選取積分步長的自動選取選定數值積分公式后，如何確定步長選定數值積分公式后，如何確定步長h以滿足給定的誤差以滿足給定的誤差 )()(122afbfhTIn梯形公式)(412nnTITInnTT2

32、用二分法只要用二分法只要其中其中fk+1/2是原分點是原分點xk,xk+1的中點的中點（記記xk+1/2）的函數值的函數值1021222nkknnfhTT且且T T2n2n可在可在T Tn n基礎上計算基礎上計算)(3122nnnTTTInTI2)2(2nnhh高斯高斯(Gauss)(Gauss)求積公式求積公式矩形公式矩形公式(1)、(2)梯形公式梯形公式(3)辛普森公式辛普森公式(4)A Ak k是與是與f f無關的常數無關的常數代數代數精度精度設設,)(kxxf用用(7)計算計算,)(badxxfI若對于若對于mk,1 ,0都有都有, IIn而當而當, 1IImkn則則稱稱In的代數精度

33、為的代數精度為m.)7()(1nkkknxfAINewton-Cotes方法方法梯形公式的代數精度（考察梯形公式的代數精度（考察T1）k=1f(x)=x222abxdxIba2)()()(21baabbfafhT3332abdxxIba2)(221baabTk=2f(x)=x2IT1IT 1梯形公式的代數精度為梯形公式的代數精度為1辛普森公式的代數精度為辛普森公式的代數精度為3高斯公式的思路高斯公式的思路取消對節點的限制，按照代數精度最大取消對節點的限制，按照代數精度最大的原則，同時確定節點的原則，同時確定節點xk和系數和系數Ak構造求積公式構造求積公式)()(22112xfAxfAG對于對于

34、11)(dxxfI使使G G2 2的代數精度為的代數精度為3 332, 1)(xxxxf)()()(221111xfAxfAdxxf確定確定2121,AAxx03/202322311222211221121xAxAxAxAxAxAAA將將f(x)f(x)代入計算得代入計算得1,3/1,3/12121AAxx)3/1 ()3/1(2ffG用用n個節點，個節點，Gn的代數精度可達的代數精度可達2n-1, 但是需解但是需解復雜的非線性方程組，實用價值不大。復雜的非線性方程組，實用價值不大。常 用 的 高 斯 公 式將將( (a,ba,b) )分小，把小區間變換為分小，把小區間變換為(-1(-1，1)

35、, 1), 再用再用G G2 2mkkkbazfzfhdxxf1)2()1 ()()(2)(322,3221) 2(1) 1 (hxxzhxxzkkkkkkmkkhaxmabhk, 1, 0,/ )(代數精度為代數精度為3節點加密時，原計算信息無法利用節點加密時，原計算信息無法利用思路思路：將積分區間分小，在小區間上用：將積分區間分小，在小區間上用n不太不太 大大的的 。而在節點加密一倍時能夠利用原節點的函。而在節點加密一倍時能夠利用原節點的函數值，可以把區間的端點作為固定節點。數值，可以把區間的端點作為固定節點。改進的高斯公式nG)()()(121bfAxfAafAGnnkkknGauss-

36、Lobatto求積公式求積公式 其中a, b為小區間的端點，nnAAxx,112為2n-2個參數，代數精度可達到代數精度可達到2n-3注意：實際計算中一般采用自適應方法確定步長注意：實際計算中一般采用自適應方法確定步長用用MATLAB 作數值積分作數值積分10nkknfhLnkknfhR1矩形矩形公式公式sum(x)輸入數組x(即fk),輸出x的和(數)cumsum(x)輸入數組x,輸出x的依次累加和(數組)梯形梯形公式公式)(2011nnkknffhfhTtrapz(x)輸入數組x,輸出按梯形公式x的積分(單位步長)trapz(x,y)輸入同長度數組 x,y,輸出按梯形公式y對x的積分(步長

37、不一定相等)用用MATLAB 作數值積分作數值積分mabhffffhSmkkmkkmn2),24(3112101220辛普森公式辛普森公式quad(fun,a,b,tol,trace)I,fn=quad()用自適應辛普森公式計算tol為絕對誤差，缺省時為10-6Gauss-Lobatto公式公式)()()(121bfAxfAafAGnnkkknquadl(fun,a,b,tol,trace)I,fn=quadl()用自適應Gauss-Lobatto公式計算 tol為絕對誤差，缺省時為10-6注意：注意：fun.m中應以自變量為矩陣的形式輸入中應以自變量為矩陣的形式輸入(點運算點運算)矩形域上計

38、算二重積分的命令：矩形域上計算二重積分的命令：dblquad(fun,xmin,xmax,ymin,ymax,tol)廣義積分、二重和三重積分長方體上計算三重積分的命令：長方體上計算三重積分的命令：triplequad(fun,xmin,xmax,ymin,ymax, zmin,zmax,tol)注：注：fun是被積函數，本身可以有自己的參數是被積函數，本身可以有自己的參數廣義積分：廣義積分：通過分析和控制誤差，轉換成普通積分通過分析和控制誤差，轉換成普通積分quadv(fun,a,b,tol,trace)向量值積分：向量值積分：用用MATLAB 作數值積分作數值積分例例. 計算計算4011s i nd xx1 1）矩形公式和梯形公式）矩形公式和梯形公式將將(0, (0, /4)100 /4)100等分等分2 2）辛普森公式和辛普森公式和Gauss-Lobatto公式公式精確、方便精確、方便無法計算用數值給出的函數的積分無法計算用數值給出的函數的積分Jifen1a.mJifen1a.mJifen1b.mJifen1b.m精確值為精確值為2數值積分的應用數值積分的應用實例實例人造衛星軌道長度人造衛星軌道長度)20(sin,costtbytax決定由短半軸長半軸rssba,21dttbtadtyxL2022222022cossin44軌道長度軌道長度yxo 近地點s1=439km