MATLAB教程：第八章 積分與數(shù)值積分

1、求符號積分函數(shù)：求符號積分函數(shù)：int格式：格式：int(f,x,a,b)功能：計算定積分功能：計算定積分dxxfba )(格式：格式：int(f,x)功能：計算不定積分功能：計算不定積分 dxxf)(使用使用int函數(shù)之前，先用函數(shù)之前，先用syms聲明聲明x是符號變量是符號變量例例8.1 求積分求積分 dxx411)1(dxxxex 102)(1)2( 32)3(2xxdx syms x y1=1/(1+x4); y2=(x*exp(x)/(1+x)2; y3=1/(x2+2*x+3); fy1=int(y1) fy2=int(y2,0,1) fy3=int(y3,-inf,+inf)二、

2、數(shù)值積分二、數(shù)值積分 在科學(xué)研究和工程技術(shù)中，經(jīng)常遇到積分的計算，在科學(xué)研究和工程技術(shù)中，經(jīng)常遇到積分的計算，雖然有些函數(shù)的不定積分可以求出其初等函數(shù)表示式，雖然有些函數(shù)的不定積分可以求出其初等函數(shù)表示式，但有更多的函數(shù)，它們的不定積分不是初等函數(shù)，但有更多的函數(shù)，它們的不定積分不是初等函數(shù)，這樣就無法利用牛頓萊布尼茲公式求出其定積分，如這樣就無法利用牛頓萊布尼茲公式求出其定積分，如,sin102 dxx,sin10 dxxx,sin2102 dtt axdxe02甚至經(jīng)常遇到只知道函數(shù)在一些離散點的值，但函甚至經(jīng)常遇到只知道函數(shù)在一些離散點的值，但函數(shù)表達式未知的情況，在上述情況下就必須以數(shù)

3、值數(shù)表達式未知的情況，在上述情況下就必須以數(shù)值方法求定積分的近似值。方法求定積分的近似值。用數(shù)值方法求定積分的近似值，通常稱為數(shù)值積分。用數(shù)值方法求定積分的近似值，通常稱為數(shù)值積分。常用的數(shù)值積分算法有：常用的數(shù)值積分算法有： 牛頓牛頓-柯特斯積分及其復(fù)化格式，自適用積分公式，柯特斯積分及其復(fù)化格式，自適用積分公式，龍貝格求積公式，高斯求積公式。龍貝格求積公式，高斯求積公式。Matlab中提供了求數(shù)值積分的函數(shù)：中提供了求數(shù)值積分的函數(shù)：(1) quad函數(shù)：基于變步長辛普森法計算積分。函數(shù)：基于變步長辛普森法計算積分。該函數(shù)調(diào)用格式：該函數(shù)調(diào)用格式：I,n=quad(fname,a,b,To

4、l,trace)其中：其中：fname是被積函數(shù)名是被積函數(shù)名 a,b是積分上下限是積分上下限 Tol是精度控制值，省卻時取是精度控制值，省卻時取0.001 Trace:控制是否顯示展現(xiàn)積分過程，取控制是否顯示展現(xiàn)積分過程，取0不展現(xiàn)不展現(xiàn) I：積分值：積分值 n：被積函數(shù)調(diào)用次數(shù)：被積函數(shù)調(diào)用次數(shù)dxxx 24sin如：求積分如：求積分ac=(x)sin(x)./xs=quad(ac,pi/4,pi/2)ac = (x)sin(x)./xs = 0.6118(2) trapz函數(shù)：函數(shù)：用梯形法計算積分，用梯形法計算積分，適用于被積函數(shù)為離散數(shù)據(jù)時，適用于被積函數(shù)為離散數(shù)據(jù)時，求函數(shù)的定積分

5、。求函數(shù)的定積分。該函數(shù)調(diào)用格式：該函數(shù)調(diào)用格式：I=trapz(x,y)dxxx 24sin如：求積分如：求積分clc,clearformat longac=(x)sin(x)./xs1=quad(ac,pi/4,pi/2)x1=pi/4:pi/50:pi/2;y1=ac(x1);s2=trapz(x1,y1)x2=pi/4:pi/100:pi/2;y2=ac(x2);s3=trapz(x2,y2)ac = (x)sin(x)./xs1 = 0.611786287055481s2 = 0.591535799623818s3 = 0.611773186493237例例8.2 神舟六號載人飛船與

6、神舟六號載人飛船與2005年年10月月12日發(fā)射升空后，在太日發(fā)射升空后，在太空中運行空中運行115小時小時32分，繞地球飛行分，繞地球飛行76圈，其先在軌道傾角圈，其先在軌道傾角42.4度、近地點高度度、近地點高度200公里、遠地點高度公里、遠地點高度347公里的橢圓軌道公里的橢圓軌道上運行上運行5圈，實施變軌后，進入圈，實施變軌后，進入343公里圓形軌道。試計算神州公里圓形軌道。試計算神州六號載人飛船在軌飛行公里數(shù)。六號載人飛船在軌飛行公里數(shù)。解解 建立坐標(biāo)系，使坐標(biāo)原點位于橢圓建立坐標(biāo)系，使坐標(biāo)原點位于橢圓中心，中心，x軸位于橢圓長軸上。軸位于橢圓長軸上。 s1=200, s2=347,

7、 地球半徑地球半徑r=6371 橢圓長半軸橢圓長半軸 arss2221 )2(2121rssa 2s1s橢圓半焦距為橢圓半焦距為c： asrc 1)(2112ssc 橢圓短半軸：橢圓短半軸： 22cab 2s1sc橢圓方程：橢圓方程： tbytaxsincosrssc )(2121)(1rs 已知已知 s1=200, s2=347, r=6371 )2(2121rssa s1=200, s2=347, r=6371 )(2112ssc 22cab 橢圓方程：橢圓方程： tbytaxsincos)2(2121rssa 1sc橢圓周長：橢圓周長： dtyxLtt 2022dttbta 202222

8、cossin4dttca 20222cos4dttaca 20222cos14 s1=200, s2=347, r=6371 )(2112ssc 22cab )2(2121rssa dttacaL 20222cos14 s1=200;s2=347; r=6371;r2=343; a=(2*r+s1+s2)/2; c=(s2-s1)/2; k=c2/a2; g=(t)sqrt(1-k*cos(t).2); L=quad(g,0,pi/2) Lch=4*a*L*5+2*pi*(r+r2)*71得：得：L1 =1.5707 Lch =3.2039e+006例例8.3 (面積的估算面積的估算)，測得某

9、人工湖邊界上，測得某人工湖邊界上n個點的坐個點的坐標(biāo)（見表），估算出該人工湖的面積標(biāo)（見表），估算出該人工湖的面積表：表：邊界上點實測坐標(biāo)邊界上點實測坐標(biāo)X2 215.015.028.628.652.052.061.861.875.475.492.992.9105.9105.9Y138.518.612.45.810.623.836.838.8Y238.578.586.881.598.6114.2160.6172.8X121.1121.1 140.8140.8 157.2157.2 180.5180.5 216.8216.8 248.5248.5 279.8279.8 310.4310.4Y12

10、5.428.718.838.650.520.45.432.8Y2176.6165.5158.6155.6186.8252.4268.8276.7X335.0335.0 360.6360.6 378.5378.5 400.1400.1 420.5420.5 448.5448.5 470.5470.5 493.1493.1Y151.956.872.866.655.565.590.5170.4Y2275.4258.2235.2228.9198.5195.5188.1170.4作出作出人工湖的平面圖：人工湖的平面圖：clc,clear,load hqdatax=A(1,:),A(4,:),A(7,:)

11、;y1=A(2,:),A(5,:),A(8,:);y2=A(3,:),A(6,:),A(9,:);fill(x,y1,g)hold onfill(x,y2,g)hold off050100150200250300350400450500050100150200250300)(2xyy )(1xyy dtxyxyA 1 .493212)()(（1 1）利用梯形法求面積）利用梯形法求面積 S=trapz(x,y2-y1); Mianji1=vpa(S,6)得得：Mianji1 = 73330.8（2）利用樣條函數(shù)積分法）利用樣條函數(shù)積分法 y=y2-y1; sc=spline(x,y); js=f

12、nint(sc); %樣條函數(shù)求積分樣條函數(shù)求積分 S=ppval(js,2,493.1); mianji2=vpa(S(2),6)得得：Mianji2 = 73586.0例例 8.4 在某山區(qū)的一個矩形區(qū)域內(nèi)，測得一些地面點的高程在某山區(qū)的一個矩形區(qū)域內(nèi)，測得一些地面點的高程(見見下表下表)，假設(shè)所觀測的矩形區(qū)域內(nèi)均是陸地，地表可近似看做一，假設(shè)所觀測的矩形區(qū)域內(nèi)均是陸地，地表可近似看做一張光滑曲面，張光滑曲面，A(0,300,420.2),B(1200,300,332.5),C(800,1000,407.5)為地表上三個點，現(xiàn)需在為地表上三個點，現(xiàn)需在A與與B，A與與C點之間的直線走向上沿

13、地點之間的直線走向上沿地表鋪設(shè)通訊線路，估計該通信線路的長度。表鋪設(shè)通訊線路，估計該通信線路的長度。 01002003004005006007008009001000110012000381.5397.3401.5396.1379.2361.5339.6320.7307.8299.6293.5285.5273.0100398.0418.5423.2416401380.5356.3336.4325.3315.5310.4306.3298.5200411.2428.3436.5433.5420.5398.5374.5353.4338.3330.2325.4323.1317.0300420.2439

14、.6448.8443.1435.5415.2392.3370.2354.1346.5342.3338.2332.5400423.8444.6456.7459.2447.6429.2407.3385.5369.4360.3356.5352.5345.2500422.1441.5459.5463.2454.6439.5419.5400.5384.5374.5369.5364.5356.4600416.2439.2454.5466.5457.6446.5429.5412.4398.5386.3381.2375.5362.1700406.3429.5442.3454.6450.2448.1435.24

15、22.2409.3399.5391.2380.2370.5800391.5414.6431.4443.9447.8440.2436.5426.5416.5406.3397.4382.2363.5900373.5395.6413.5426.6433.5434.5430.5424.5416.5407.6390.0388.6370.01000352.3374.5391.6405.5414.3418.5416.8413.6407.5398.5389.5376.5350.5區(qū)域區(qū)域0,1200;0,1000內(nèi)觀測點處的高程表內(nèi)觀測點處的高程表 （單位：米）（單位：米） 記地表曲面方程為記地表曲面方程為z

16、=f(x,y),由于由于A與與B的直線走向平行于的直線走向平行于x軸，軸，沿該走向的地表曲線沿該走向的地表曲線LAB為：為： 300)300,(yxfz曲線曲線LAB的弧長為：的弧長為：dxgx 120002)(1)(xg 問題求解問題求解: (1)通過插值法求出函數(shù)通過插值法求出函數(shù)z=g(x) (2)通過數(shù)值積分求出該地表曲線的長度通過數(shù)值積分求出該地表曲線的長度地表高程數(shù)據(jù)存于文件：地表高程數(shù)據(jù)存于文件：data804clear,clcload data804x=0:100:1200;y=0:100:1000;g=z(4,:);pp=csape(x,g);dp=fnder(pp);gx=

17、ppval(dp,x);hwf=sqrt(1+gx.2);s=trapz(x,hwf);dAB=vpa(s,8)得得：dAB =1211.9855 300)300,(yxfzdxgx 120002)(1)(xg A與與B點的連線點的連線(y=300)的的數(shù)據(jù)位于地面高程矩陣的數(shù)據(jù)位于地面高程矩陣的第四行上，第四行上，記地表曲面方程為記地表曲面方程為z=f(x,y),A(0,300,420.2)與與C(800,1000,407.5)的直線走向與的直線走向與x軸正軸正向夾角為向夾角為u，則：，則：08003001000tan u,87 ,1138cos u,1137sin u記水平面上記水平面上A

18、1(0,300)與與C1(800,1000)連線上任一點連線上任一點P(x,y)到點到點A1的距離為的距離為t，有，有,cosutx utysin300 沿該走向的地表曲線沿該走向的地表曲線LAC為：為： cos sin300)sin300,cos(utxutyututfz)(th )1131000( t沿該走向的地表曲線沿該走向的地表曲線LAC為：為： cos sin300)sin300,cos(utxutyututfz)1131000( t,1138cos u,1137sin u曲線曲線LAC的弧長為：的弧長為：dthyxttt 1130222)()()()(th dtht 11302)(1load data804x=0:100:1200;y=0:100:1000;t=0:100*sqrt(113);xx=8/sqrt(113)*t;yy=300+7/sqrt(113)*t;x0,y0=meshgrid(x,y);h=interp2(x0,y0,z,xx,yy,spline);pp1=csape(t,h); %樣條插值求函數(shù)樣條插值求函數(shù)h(t)dp1=fnder(pp1); %求樣條插值函數(shù)求樣條插值函數(shù)h(t)的導(dǎo)數(shù)的導(dǎo)數(shù)ht=fnval(dp1,t);hwf=sqrt(1+ht.2); %求出被積函數(shù)求出被積函