数值积分的matlab实现5551

.实验10数值积分实验目的：1．了解数值积分的基本原理；2．熟练掌握数值积分的MATLAB实现；3．会用数值积分方法解决一些实际问题。实验内容：积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所以仍然得不到积分的精确值，如不定积分1sinxdx。这时我们一般考虑用数值方法计算其x0近似值，称为数值积分。10.1数值微分简介yf(x)在设函数x*可导，则其导数为f(x*)limf(x*h)f(x)*（10.1）hh0yf(x)以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得如果函数其近似值f(x*)f(xh)f(x)**（10.2）h表10-1xxyxx…………xy0122ny0yy1nyf(x)在一般的，步长越小，所得结果越精确。（10.2）式右端项的分子称为函数hx*的差分，分母称为自变量在x*的差分，所以右端项又称为差商。数值微分即用差商近似.代替微商。常用的差商公式为：f(xh)f(xh)f(x)（10.3）（10.4）002h3y4yy20f(x)0012hf(x)y4y3yn（10.5）n2n12hn其误差均为O(h2)，称为统称三点公式。10.2数值微分的MATLAB实现MATLAB提供了一个指令求解一阶向前差分，其使用格式为：dx=diff(x)xx,xx,,xx,这样基于两点的数值导其中x是维数组，dx为1维数组nn2132n1数可通过指令diff(x)/h实现。对于三点公式，读者可参考例1的M函数文件diff3.m。例1用三点公式计算yf(x)在x1.0,1.2,1.4处的导数值，f(x)的值由下表给出。x1.01.11.21.31.4f(x)0.25000.22680.20660.18900.1736解：建立三点公式的M函数文件diff3.m如下：functionf=diff3(x,y)n=length(x);h=x(2)-x(1);f(1)=(-3*y(1)+4*y(2)-y(3))/(2*h);forj=2:n-1f(j)=(y(j+1)-y(j-1))/(2*h);endf(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*h);.在MATLAB指令窗中输入指令：x=[1.0,1.1,1.2,1.3,1.4];y=[0.2500,0.2268,0.2066,0.1890,0.1736];diff3(x,y)运行得各点的导数值为：-0.2470，-0.2170，-0.1890，-0.1650，-0.0014。所以x1.0,1.2,1.4处的导数值分别为-0.2470，-0.1890和-0.0014。yf(x)在对于高阶导数，MATLAB提供了几个指令借助于样条函数进行求导，详细使用步骤如下：step1：对给定数据点（x,y），利用指令pp=spline(x,y)，获得三次样条函数数据pp，供后面ppval等指令使用。其中，pp是一个分段多项式所对应的行向量，它包含此多项式的阶数、段数、节点的横坐标值和各段多项式的系数。step2：对于上面所求的数据向量pp，利用指令[breaks,coefs,m,n]=unmkpp(pp)进行处理，生成几个有序的分段多项式pp。step3：对各个分段多项式pp的系数，利用函数ppval生成其相应导数分段多项式的系数，再利用指令mkpp生成相应的导数分段多项式step4：将待求点xx代入此导数多项式，即得样条导数值。上述过程可建立M函数文件ppd.m实现如下：functiondy=ppd(pp)[breaks,coefs,m]=unmkpp(pp);fori=1:mcoefsm(i,:)=polyder(coefs(i,:));enddy=mkpp(breaks,coefsm);于是，如果已知节点处的值x,y，可用下面指令计算xx处的导数dyy：

.pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx);例2基于正弦函数ysinx的数据点，利用三点公式和三次样条插值分别求导，并与解析所求得的导数进行比较。解:编写M脚本文件bijiao.m如下：h=0.1*pi;x=0:h:2*pi;y=sin(x);dy1=diff3(x,y);pp=spline(x,y);dy=ppd(pp);dy2=ppval(dy,x);z=cos(x);error1=norm(dy1-z),error2=norm(dy2-z)plot(x,dy1,'k:',x,dy2,'r--',x,z,'b')运行得结果为：error1=0.0666，error2=0.0025，生成图形见图10.1。图10.1三点公式、三次样条插值与解析求导比较图显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图2.15可知利用三次样条插值求导所得曲线与解析求导曲线基本重合，而三点公式在极值点附近和两个端点附近误差较大，其它点吻合的较好。10.3应用示例：湖水温度变化问题问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低。这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。.对某个湖的水温进行观测得数据见表10-2。表10-2某湖的水温观测数据深度（m）温度（℃）02.34.99.113.713.918.311.722.911.127.211.122.822.822.820.6试找出湖水温度变化最大的深度。1．问题的分析湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。对于给定的数据，可以利用数值微分计算各深度的温度变化值，从而得到温度变化最大的深度，但考虑到所给的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，以得到更准确的结果。2．模型的建立及求解hTTT(h)，并假定函数T(h)可记湖水的深度为（m），相应的温度为（℃），且有导。对给定的数据进行三次样条插值，并对其求导，得到T(h)的插值导函数；然后将给定的深度数据加密，搜索加密数据的导数值的绝对值，找出其最大值及其相应的深度，相应的MATLAB指令如下：h=[02.34.99.113.718.322.927.2];T=[22.822.822.820.613.911.711.111.1];hh=0:0.1:27.2;pp=spline(h,T);dT=ppd(pp);dTT=ppval(dT,hh);[dTTmax,i]=max(abs(dTT)),hh(i)plot(hh,dTT,'b',hh(i),dTT(i),'r.'),gridon运行得导函数绝对值的最大值点为：=11.4，最大值为1.6139,即湖水在深度为11.4mh.时温度变化最大，如图10.2所示（黑点为温度变化最大的点）。图10.2湖水温度变化曲线图10.4数值积分简介考虑定积分bf(x)dx（10.6）afx如果被积函数()是以列表形式给出，则其求解思想同数值微分类似，即用逼近多项式()近似地代替被积函数()，然后计算积分PxnfxbP(x)dx，得（10.6）式的近似值；na如果被积函数的原函数不是初等函数，则将积分区间进行细分，对每个小区间，用一个近似fx函数代替被积函数()，然后积分得（10.6）式的近似值。这两种类型最终都可归结为函fx数()在节点xk上的函数值f(x)的某种线性组合，即下面数值求积公式：kIf(x)dx或（10.7）nAf(x)bkkak0Af(x)Rf（10.8）bf(x)dxnkkak0其中Rf为截断误差。此误差可用代数精度衡量，代数精度越高，误差越小；反之误差越大。代数精度是用来衡量数值积分公式近似程度的办法，如果()是一个次数不超过的fxmfx代数多项式，（10.7）式等号成立；而当()是一个m1次多项式时，（10.7）式不能精成确立，则称（10.7）式的代数精度为m。.选取不同的近似函数，可产生不同的数值求积公式，常见的有：梯形公式、辛普森公式和高斯公式。10.5数值积分的MATLAB实现MATLAB提供了下面几个函数计算积分，其使用格式分别为：采用梯形公式计算积分（h1）,x为f(k0,1,,)n（1）trapz(x)k（2）quad('fun',a,b,tol)（3）quadl('fun',a,b,tol)采用自适应Simpson法计算积分采用自适应Gauss-Lobatto法计算积分其中fun为被积函数；tol是可选项，表示绝对误差，a，b为积分的上、下限。例1分别利用梯形公式、Simpson公式和Gauss-Lobatto法计算112dx，并与x0其精确值比较。解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方法求得的数值解比较，其MATLAB指令为：symsxxz0=simple(int('sqrt(1+xx^2)',0,1))z=double(z0);z=vpa(z,8)x=0:0.01:1;y=sqrt(1+x.^2);z1=trapz(y)*0.01;z1=vpa(z1,8),err1=z-z1;err1=vpa(err1,8)z2=quad('sqrt(1+x.^2)',0,1);z2=vpa(z2,8),err2=z-z2;err2=vpa(err2,8)z3=quadl('sqrt(1+x.^2)',0,1);z3=vpa(z3,8),err3=z-z3;err3=vpa(err3,8)1(2ln(21))1.1477936，三种公式计算得数值积分值分别为运行得精确值为21.1477995，1.1477935和1.1477936，其相应误差分别为-.59e-5，.1e-6和0.，由三者误差可见，Gauss-Lobatto法计算最为精确，Simpson公式次之，梯形公式最差，但它也.能精确到小数点后5位数。例2人造地球卫星轨道可视为平面上的椭圆。我国第一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为6371km，求该卫星的轨道长度。t2),a,b分别是长、短半xacost,ybsint(0解：卫星轨道椭圆的参数方程为轴，则根据所给数据知=6371+2384=8755,=6371+439=6810。ab由对弧长的曲线积分知识知，椭圆的长度为1L4(a2sin2tb2cos2t)2dt20上积分称为椭圆积分，它无法用解析方法计算，可用计算其数值解，编写M函数文件如下：functiony=y(t)a=8755;b=6810;y=4*sqrt(a^2*sin(t).^2+b^2*cos(t).^2);在MATLAB指令窗中输入以下指令：l=quad('y',0,pi/2)运行得结果为：l=4.9090e+004。即卫星轨道长度为49090km。对于用列表形式给出的函数，上述方法不再适用，可利用指令spline构造三次样条插值函数，再计算积分，具体步骤可参考例2。例3在桥梁的一端每隔一段时间记录1min有几辆车过桥，得到数据见表10-3：表10-3过桥车辆数据0:0022:0024:0005:0026:0057:0088:0025车辆数/辆车辆数/辆车辆数/辆9:001210:30511:301012:301214:00716:00917:002818:002219:001020:00921:001122:00823:00924:003试估计一天通过桥梁的车流量。.解：记记录时刻为时，相应的车辆数为x(t)，则所求车流量即为计算积分24tx(t)dt，0则在MATLAB指令窗中输入下面指令：x=[0,2,4,5,6,7,8,9,10.3,11.3,12.3,14,16,17,18,19,20,21,22,23,24];y=[2,2,0,2,5,8,25,12,5,10,12,7,9,28,22,10,9,11,8,9,3];pp=spline(x,y);s1=quadl('fun',0,24,[],[],pp)%利用三次样条插值计算积分s2=trapz(x,y)%利用梯形公式计算积分其中M函数文件fun.m为：functionvf=fun(x,pp)vf=ppval(pp,x);运行得三次样条插值计算所得车流量为212辆，梯形公式计算所得车流量为216辆。10.6数值积分的建模示例：煤炭储量计算问题问题：某煤矿为估计其煤炭的储量，在该矿区内进行勘探，得到数据如表10-4。试估算出该矿区（1xy4,15）煤炭的储量。表10-4某煤矿勘探数据表编号12345678910x坐标（km）11121314152122232425y坐标（km）h煤炭厚度（m）13.7225.808.471325.2722.3215.4721.3314.4924.8326.19编号111214151617181920x坐标（km）31323334354142434445y坐标（km）h煤炭厚度（m）23.2826.4829.1412.0414.5819.9523.7315.3518.0116.291．问题的分析问题给出了很多点对应的煤炭厚度，显然整个煤矿可以看作是一个巨大的曲顶柱体（见图10.3，此图经过插值得到），而煤炭的储量即为此立体的体积。要计算此立体的体积，可.以利用插值得到曲顶柱体的顶面函数，再对其积分；也可以将此曲顶柱体分割成若干个细的曲顶柱体，用数值方法计算这些细曲顶柱体的体积，再对其求和即得原曲顶柱体的体积。图10.3煤炭厚度曲面图2．模型的建立及求解以煤炭的厚度为三维空间中的坐标建立空间坐标系。记煤炭的厚度为，则它是坐标zhx,y的二元函数，即zW(x,y)，则由二重积分的知识可知，此煤矿的煤炭储量为W(x,y)dxdy（10.9）D其中D{(x,y)|1x4,1y5}。由于函数(x,y)只给出了一些离散点上的函数值，无法直接计算上述二重积分，所以下面采用数值积分的方法计算其值。由数值积分的知识知，计算定积分有复合梯形公式为bf(x)dxh[f(a)f(b)2n1f(x)]k（10.10）（10.11）2ak1其中为步长，hx(k0,1,,)n为节点，且有xakh。kk由（10.9）式得bg(x)dxdWxyyx(,)ddbaca.其中g(x)(x,y)dy，则由（10.10）式可得dchx[g(a)g(b)2Wbg(x)dxg(x)]jn1（10.12）2aj1而(a,y)dyy[(a,c)(a,d)2(a,y)]hm1g(a)g(b)d2kck1h(b,y)dyy[(b,c)(b,d)2m1(b,y)]d2kck1hg(x)k(x,y)dyy[(x,c)(x,d)2m1(x,y)]d2kkkkkck1所以有Whhm1(a,y)(b,y)[(a,c)(a,d)(b,c)(b,d)2xy4kkk12[(x,c)(x,d)2m1(x,y)]]n1jjjk（10.13）j1k1hh[5(a,c)5(a,d)5(b,c)5(b,d)xy42m(a,y)(b,y)2n[(x,c)(x,d)]4nm(x,y)kkjjjkk0j0j0k0考虑到给定的数据较少，由此产生的误差较大，所以利用插值后的数据计算（10.13）式，相应的MATLAB计算指令如下：x=[11111222223333344444]*1000;y=[12345123451234512345]*1000;z=-[13.7225.808.4725.2722.3215.4721.3314.4924.8326.1923.2826.4829.1412.0414.5819.9523.7315.3518.0116.29];hx=10;hy=10;cx=1000:hx:4000;cy=1000:hy:5000;[X,Y]=meshgrid(cx,cy);n=length(cx);m=length(cy);.Z=griddata(x,y,z,X,Y,'v4');surf(X,Y,Z)%插值%绘制图10.3W=-hx*hy*(-5*Z(1,1)-5*Z(1,n)-5*Z(m,1)-5*Z(m,n)+2*(sum(Z(1,:)+Z(m,:))+sum(Z(:,1)+Z(:,