2019年数值积分的matlab实现

1、f(x9h)?f(x?h)?oo?(x)f o2h?3y?4y?y?2io?)(fx y?4y?3y?Xmm“?x)f(三点公式。其误差均为102数值微分的MATLAB实现实验10数值积分实验目的：1. 了解数值积分的基本原理：2. 熟练掌握数值积分的MATLAB实现：3. 会用数值积分方法解决-些实际问题。实验内容：积分是数学中的个基本概念，在实际问题中也有很广泛的应用。同微分样，在微积分中， 它也是通过极限定义的，由于实际问题中遇到的函数般都以列农形式给岀，所以常常不能用来 宜接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所SlnX 1?Xd这时我们般考虑用数值方法计算其如

2、不定积分。以仍然得不到积分的精确值， Xo近似值，称为数值积分。10. 1数值微分简介x)f(x?y在设函数可导，则其导数为-)x?f(x(?h)f*?lim?(Xf) )(10. 1 homy?f(x)以列农形式给出(见表10-1),则其精确值无法求得，但可由 下式求得如果函数其近似值)(fk?h)f(x?xf)( 10. 2)(h 表 10-1Xyhy?f(x)在(10.2)式右端项的分了称为函数般的，步长越小，所得结果越精确C= xx的差分， 所以右端项又称为差商。的差分，分母称为自变量在数值微分即用差商近似代替微商。常用的差 商公式为：(10.3 ) (10.4)oh2.(10.5)

3、n21)h(,称为统称MATLAB提供了 个指令求解阶向前差分，其使用格式为: dx=diff(x)?nx?,XXX,x?X?l?n,这样展于两点的数值导dx维数组，为其中X是维数组】述1数可通过指令diff(x)h实现。对于三点公式，读者可参考例1的M函数文件diff3.mox7f(xx)yf的值由下衣给用三点公式计算处的导数值，在1. 0,1. 2,1. 41例X)(xf1.01.31. 40. 25000. 226S0. 20660. 1S900.1736解:建立三点公式的M函数文件diff3.m如下:Rinction f=df(x,y)n=length(x) ;h=x(2)-x( 1)

4、;f(lH-3y(l)÷4*y(2)-y(3)(2*h);forj=2:n-lf(j)=(y(ji)-y(j-i)y(2*h);endf(nXy(n-2)-4y(n4)÷3*y(n)(2*h);在MATLAB指令窗中输入指令：x=1Q1.1,1.2,L3,14; y=0.2500,0.2268,0.2066,0.1890,0.1736;diff3(%y)y?f(x)所以，-0.0014。，运行得各点的导数值为：-0.2470, -0.2170-0. 1890,-0. 1650x71. 0, 1. 2, 1. 4 处的导数值分别为-0 2470> -0. 1890 和-

5、0. 0014 在。对于高阶导数，MATLAB提供了几个指令借助于样条函数进行求导，详细使用步骤如下：stepl：对给定数据点(x, y),利用指令PP=SPIine(x, y),获得三次样条函数数据pp,供后面 PPVaI等指令使用。其中，PP是个分段多项式所对应的行向量，它包含此多项式的阶数、段数、 节点的横坐标值和各段多项式的系数。step2:对于上而所求的数据向 PP»利用指令breaks, coefs, m, n=unmkpp(pp)进行处理，生 成几个有序的分段多项式PP。step3:对各个分段多项式PP的系数，利用函数PPVaI生成其相应导数分段多项式的系数，再利 用指

6、令mkpp生成相应的导数分段多项式step4:将待求点代入此导数多项式，即得样条导数值。上述过程可建立M函数文件ppd. m实现如下：ficton dy=d(p)breaks.coefs,m=mmk();for=l :mCOefSm(ij)=olyder(coefs(,:);enddy=nk(breaks,coefsm);于是，如果已知节点处的值x,y,可用下面指令计算XX处的导数dyy：PP=SPIine(X,y).dr=ppd(pp):dy3r=ppval(dy.xx);y?SlnX的数据点，利用三点公式和三次样条插值分别求导，并与基于正弦函数例2解析所求得 的导数进行比较。解:编写M脚本

7、文件bi jiao, m如下：h=0.1*i=0*iir=sm(x); dyl=d3(x,y); PP=SPlme(X,y)：dy=ppd(pp):dy2=ppval(dy.x); Z=COS(x);error 1 =norm(dy 1 -z).error2=norm(dy2-z) )rb,r-,.z.lot(x,dyl,k x,dy2,运行得结果为:errorl =0. 0666, error2 =0. 0025,生成图形见图10IO0.5 -0 -0.5 -1o三点公式、三次样条插值与解析求詁比较图图10.1显然利用三次样条插值求导所得谋差比三点公式求导小很多，同时由图2. 15可知利用三

8、次样条 插值求导所得曲线与解析求导曲线基本重合,而三点公式在极值点附近和两个端点附近谋差较人, 其它点吻合的较好。10.3应用示例：湖水温度变化问题问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低。 这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。对某个湖的水温 进行观测得数据见表10-2 o表10-2某湖的水温观测数据深度()02.34.99.113.71S. 322.927.2温度(C)22. S22. S22. S20.613.911.711. 111.1试找出湖水温度变化最人的深度。1. 问题的分析湖水的温度可视为关于深度的函数，

9、于是湖水温度的变化问题便转化为温度函数的导数问题，显 然导函数的最人绝对值所对应的深度即为温度变化最人的深度。对于给定的数据，但考虑到所给 从而得到温度变化最大的深度，可以利用数值微分计算各深度的温度变化值，的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，以得 到更准确的结果。2. 模型的建立及求解TT7T(h)T(h)可(°C),相应的温度为，且有记湖水的深度为，并假定函数(m)导。T(h)的插值导函数：然后将给定对给定的数据进行三次样条插值，并对其求导，得到的深度数据 加密,捜索加密数据的导数值的绝对值，找出其最人值及其相应的深度，相应的MATLA

10、B指令如 下：h=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1;hh=O:0.1:27.2;PP=SPIme(h.T):dT=ppd(pp):dTT=ppval(dT.hli);dTTmax.=max(abs(dTT),hh()').grd Onplot(hh,dTThh().dTT(),b,r.h= 11. 4,最人值为运行得导函数绝对值的最人值点为：1.6139,即湖水 在深度为11.4m时温度变化最人，如图10.2所示(黑点为温度变化最人的点。0.5图10. 2湖水沮度变化曲

11、线图10.4数值积分简介考虑定积分b?f(x)dx(10. 6) df(x)是以列衣形式给出，则其求解思想同数值微分类似，即用逼近多项如果彼积函数b?P(x)dx)P(x)(xf, 得(式10.6,然后计算积分近似地代替被积函数)式的近似值；皿如果彼积函数的原函数不是 初等函数，则将积分区间进行细分，对每个小区间，用个近似f(x),然后积分得(10.6)式的 近似值。这两种类型最终都可归结为函函数代替被积函数Xf(X)f(x的某种线性组合，即下面数 值求积公式：在节点上的函数值数Iinb?)(XX(X)d?fAIf ) 10.7或(kka"knb?PRAf(x)?f(x)dx)10.

12、8(IiaO匕?fR 为截断谋差。此谋差可用代数精度衡量，代数精度越高，谋差越小：反之误其中差越人。m)(xf的是个次数不超过代数精度是用来衡量数值积分公式近似程度的办法，如果lm?)f(x)式不能精是 个(代数多项式，(10. 7)式等号成立：而当10.7次多项式时，m。确成立，则称(10. 7) 式的代数精度为选取不同的近似函数，可产生不同的数值求积公式，常见的有：梯形公式、辛普 森公式和高斯公式。MATLAB实现10. 5数值积分的MATLAB捉供了下面几个函数计算积分，其使用格式分别为：)nO,lf(klh?»?,x) 1 O trapz(X)采用梯形公式计算积分(k法计算积

13、分)quad('fun', a, b, tol)采用自适应SimPSOn (2法计算积分)QUadI (, fun', a, bf tol)釆用自适应GaUSS-LObattO (3 a, b为积分的上、下限。其中fun为被积函数：tol是可选项，农示绝对谋差,i2?x?XdL并与其法计算GaUSS-LObattO 例1分别利用梯形公式、SilnPSOn公式和o精确值比较。解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方 法求得的数值解比较，其MATLAB指令为：SymSXXZO=Suile(int(,sqrt(l+xxA2)0.1)z=doubl

14、e(z0);Z=VPQ(Z,8)X=OrOOl :1; y=sqrt(l+x 人2);zl=traz(y)*0.0Ezl=rpa(zL8).el=z-zlerrl=,a(errl,8) z2=quad('sqrt(l+xA2)0.1);z2=pa(z2,8),err2=z-z2;err2=rpa(err2.8) z3=quadlCsqft(l+x.A2y,0,l);z3=pa(z3,8Xen3=z-z3;err3=pa(err3,8)l(2?ln(2?l)?l. 1477936,三种公式计算得数值积分值分别为运行得精确值为一 21.1477995,1. 1477935和1. 14779

15、36,其相应误差分别为-59e5, . le-6和0.,由三者误差可见, GaUSS-LObattO法计算最为持确，SimPSOn公式次之，梯形公式最差，但它也能持确到小数点后 5位数。例2人造地球卫星轨道可视为平面上的椭圆。我国第颗人造地球卫星近地点距地球农面439km, 远地点距地球农面2384km,地球半径为6371km,求该卫星的轨道长度。?),a,2(0t?t?bsm?,cos?Xatyb分别是长、卫星轨道椭圆的参数方程为解:短半ab=6371+439=68IOo 轴，则根据所给数据知=6371+2384=8755,由对弧长的曲线积分知识知，椭圆的长度为?12222?tdcostL?

16、4)(asmt?b22o上积分称为椭圆积分，它无法用解析方法计算，可用计算其数值解，编写H函数文件如下：Rinction y=y(t)a=8755:b=6810;)'=4*sqrt(aA2*sin(t).A2-bA2*cos(t).A2);在MATLAB指令窗中输入以下指令：l=quad('y'.0.pi2)运行得结果为：1=4. 9090e+004即卫星轨道长度为4909OknlO对于用列农形式给出的函数，上述方法不再适用，可利用指令SPIine构造三次样条插值函数， 再计算积分，具体步骤可参考例2。例3在桥梁的端每隔段时间记录Imin有几辆车过桥，得到数据见<

17、 10-3：表10-3过桥车辆数据时间0:002:004:005:006:007:008:00辆车辆数/29025S25时间9:0010:3011:3012:3014:0016:0017:00辆车辆数/12 O10127928时间18:0019:0020:0021:0022:0023:0021:00辆车辆数/10911893试估计天通过桥梁的车流量。24?tx(t)dt)(tx,解：记记录时刻为时，相应的车辆数为，则所求车流量即为计算积分0MATLAB 则在指令窗中输入下面指令:x=02,4,5,6,7,8,9,103,113,123,14,16,17,18,19,20,21,22,2324；

18、 y=22025,8,25,12,5,10.12,7.9,2&22,10911,893;% 利用三次样条插值计算积分 PP=SPIme(X.y)isl=quadl(fun0,24,.)s2=trapz(x,y)%利用梯形公式计算枳分其中H函数文件fun.m为：fiction XFf=ftn(x.)Vrf=PPVai(PP.x);运行得三次样条插值计算所得车流量为212辆，梯形公式计算所得车流量为216辆。10.6数值积分的建模示例：煤炭储量计算问题问题:某煤矿为估计其煤炭的储量在该矿区内进行勘探,得到数据如衣10-4试估l?x",l?y?5) 煤炭的储量。算出该矿区(表IoT

19、某煤矿勘探数据表編号193467S910X坐标(km1111122222y坐标(km1234 O1234511 ) (m煤炭师度13. 7225. 80S. 4725. 2799 3°15.4721.3311. 4924. S326. 19編号11121314151617IS1920X ) km坐标(3333344441y坐标(km12341231511 (In)煤炭伸度23. 2826. 4S29. 1412. 0414. 5819.9523. 7315. 351S.0116. 291.问题的分析问题给出了很多点对应的煤炭厚度，显然整个煤矿可以看作是个巨人的曲顶柱体(见图10.3,

20、 此图经过插值得到)，而煤炭的储虽即为此立体的体积。要计算此立体的体积，可以利用插值得 到曲顶柱体的顶而函数，再对其积分；也可以将此曲顶柱体分割成若干个细的曲顶柱体，用数值 方法计算这些细曲顶柱体的体积，再对其求和即得原曲顶柱体的体积。1000 1000图10.3煤炭M度曲而图2.模型的建立及求解Zli.则它是坐标坐标建立空间坐标系。记煤炭的厚度为以煤炭的厚度为三维空间中的 ?W)y(xZEyX,为的二元函数，即，则由二重积分的知识可知，此煤矿的煤炭储量?(x，y)dxdyW?(10.9)D?|D?(x,y)l?x?4,ly?5其中。7)yx,(H给出了 些离散点上的函数值，无法直接计算上述二

21、重 积分，所以由于函数下面采用数值积分的方法计算其值。由数值积分的知识知，计算定积分有复合梯形公式为 nnhb?f(x)(b)?2df(x)xf(a)f(10. 10)k2akn?)n 1 ,x(k01ikhxa»?O为步长，其中为节点，且有吐由(10. 9)式得Mb? ? ? ? 7xdydx9g(x)d? WyX(J (10.11)?®a?ydy)(xg(x)?,)式可得其中，则由(10. IOehIQ'b?xg(x?2)(a?)?g(b)g(x)dx?gW(10. 12)j2aj?l 而hm“!d?y?,)?2y)(g(a)?(a,c)?a(a,d(a,y)d

22、y?一I叱krid?y?2y)(b,(b,c)?(b,dg(b)?)bCy)dy? _k2cik?hm?id?y? ?)(x,)?y(x,dg(x)?)(x,y)dy?2x(,cu±kkk2ci?k 所以有hhnri?Y?)yb)(a,y)?(b,c)?(b0)W?2Q(a,c?(aqMiwi加？?y?2)(xc(xj?(Xq)2?顷ij( 10. 13) hli?(b,dc)?5)(a,d?)?5(b,?5(a,c)?5 43?区丫)(人(1)勺(b,)?24(x,c)?2?y(a,)呻Jg 考虑到给定的数据较少，由此产生的误差较人，所以利用插值后的数据计算(10. 13)式，相应

23、的MATLAB计算指令 如下：x=l 111122222333334444 4n000;y=l 234512345123451234 5*1000:z=-13.72 25.80 8.47 25.27 22.32 15.47 21.33 14.49 24.83 26.19 23.28 26.48 29.14 12.04 14.5819.95 23.73 1535 18.01 16.29;hx=10:hy=10;CX=IoO(xW000;Cy=IoOOhy5000;X, Y= meShgnd(CXXy);II=Iength(Cx);m=length(cy);Z=gndd 也(xy,zXY,PV);%插值Surf(X YZ)% 绘制图 10.3W=-hx*hy*(-5*Z(l,l)-5*Z(l,n)-5*Z(md)-5*Z(m.n)+2*(suni(Z(l,:)*Z(mj)+sum(Z(:,l)+Z(:>n)+4*sum(sum(Z) ”4w：s.ml,即煤矿的煤炭储虽约为2. 5242X10运行得X=N 5242实验任务：D?D(t)如下农所示：1. 一个物体的运动距离tS.09.010.011.012.0OtD1