数值积分的matlab实现

(10.2)(10.2)....3．会用数值积分方法解决一些实际问题。以仍然得不到积分的精确值，如不定积分j1sinxdx。这时我们一般考虑用数值方法计算其0x近似值，称为数值积分。设函数y=f(x)在x*可导，则其导数为f,(x*)=limf(x*+h)f(x*)h0h .1)如果函数y=f(x)以列表形式给出(见表10-1)，则其精确值无法求得，但可由下式求得其近似值f,(x*)必f(x*+h)f(x*)h表10-1xyxxxn2xyxxxn2yyy2n2一般的，步长h越小，所得结果越精确。(10.2)式右端项的分子称为函数y=f(x)在x*的差分，分母称为自变量在x*的差分，所以右端项又称为差商。数值微分即用差商近似用的差商公式为：2hf(x)必2h 0.4)......f(x)必n2h .5)其误差均为O(h2)，称为统称三点公式。MATLAB提供了一个指令求解一阶向前差分，其使用格式为：dx=diff(x)2132n1m例1用三点公式计算y=f(x)在x=1.0,1.2,1.4处的导数值，f(x)的值由下表给xf(x)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);f(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*h);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)step1：对给定数据点(x,y)，利用指令pp=spline(x,y)，获得三次样条函数数据pp，的阶数、段数、节点的横坐标值和各段多项式的系数。step2：对于上面所求的数据向量pp，利用指令[breaks,coefs,m,n]=unmkpp(pp)进行段多项式pp。functiondy=ppd(pp)[breaks,coefs,m]=unmkpp(pp);........fori=1:mcoefsm(i,:)=polyder(coefs(i,:));dy=mkpp(breaks,coefsm);xyxx处的导数dyy：pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx);例2基于正弦函数ysinx的数据点，利用三点公式和三次样条插值分别求导，并与解析所求得的导数进行比较。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')ror显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图2.15可知利用三次样条插值求导所得曲线与解析求导曲线基本重合，而三点公式在极值点附近和两个端点附近误差较大，其它点吻合的较好。问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低。这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。个湖的水温进行观测得数据见表10-2。深度(m)温度(℃)0试找出湖水温度变化最大的深度。湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导........的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，2．模型的建立及求解记湖水的深度为h(m)，相应的温度为T(℃)，且有T=T(h)，并假定函数T(h)可对给定的数据进行三次样条插值，并对其求导，得到T(h)的插值导函数；然后将给定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时温度变化最大，如图10.2所示(黑点为温度变化最大的点)。考虑定积分jbfxdx(10.6)a如果被积函数f(x)是以列表形式给出，则其求解思想同数值微分类似，即用逼近多项式P(x)近似地代替被积函数f(x)，然后计算积分jbP(x)dx，得(10.6)式的近似值；anna函数代替被积函数f(x)，然后积分得(10.6)式的近似值。这两种类型最终都可归结为函数f(x)在节点x上的函数值f(x)的某种线性组合，即下面数值求积公式：kkI=jbf(x)dxnAf(x)或(10.7)akkak=0......aak=0 .8)其中R[f]为截断误差。此误差可用代数精度衡量，代数精度越高，误差越小；反之误代数多项式，(10.7)式等号成立；而当f(x)是一个m+1次多项式时，(10.7)式不能精确成立，则称(10.7)式的代数精度为m。选取不同的近似函数，可产生不同的数值求积公式，常见的有：梯形公式、辛普森公式MATLAB提供了下面几个函数计算积分，其使用格式分别为： (1)trapz(x)采用梯形公式计算积分(h=1),x为f(k=0,1,,n)k(2)quad('fun',a,b,tol)采用自适应Simpson法计算积分(3)quadl('fun',a,b,tol)采用自适应Gauss-Lobatto法计算积分0解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方法求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，三种公式计算得数值积分值分别为2轴，则根据所给数据知a=6371+2384=8755,b=6371+439=6810。由对弧长的曲线积分知识知，椭圆的长度为......1L=4j2(a2sin2t+b2cos2t)2dt0上积分称为椭圆积分，它无法用解析方法计算，可用计算其数值解，编写M函数文件如下：functiony=y(t)a=8755;b=6810;y=4*sqrt(a^2*sin(t).^2+b^2*cos(t).^2);l=quad('y',0,pi/2)对于用列表形式给出的函数，上述方法不再适用，可利用指令spline构造三次样条插值函数，再计算积分，具体步骤可参考例2。例3在桥梁的一端每隔一段时间记录1min有几辆车过桥，得到数据见表10-3：22025857924:00车辆数/辆9893试估计一天通过桥梁的车流量。解：记记录时刻为t时，相应的车辆数为x(t)，则所求车流量即为计算积分j24x(t)dt，0x=[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)%利用梯形公式计算积分functionvf=fun(x,pp)vf=ppval(pp,x);储量计算问题4。试估算出该矿区(1x4,1y5)煤炭的储量。123456789x坐标(km)1111122222y坐标(km)123451234526.19x坐标(km)3333344444........y坐标(km)1234512345煤炭厚度h(m)23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29图10.3，此图经过插值得到)，而煤炭的储量即为此立体的体积。要计算此立体的体积，可曲顶柱体，用数值方法计算这些细曲顶柱体的体积，再对其求和即得原曲顶柱体的体积。2．模型的建立及求解zh(10.9)(10.9)D由于函数Q(x,y)只给出了一些离散点上的函数值，无法直接计算上述二重积分，所以下面采用数值积分的方法计算其值。由数值积分的知识知，计算定积分有复合梯形公式为jbf(x)dx~h[f(a)+f(b)+2f(x)]a2k(10.10)kk由(10.9)式得(10.11)(10.11)......其中g(x)=jdQ(x,y)dy，则由(10.10)式可得c(10.12)W=jbg(x)dxhx[g(a)+g(b)+2g(x)](10.12)a2jj=1而c2kk=1c2kk=1kck2kkkkk=1有4kkjjjkj=1k=1(10.13)kkjjjkk=0j=0j=0k=0考虑到给定的数据较少，由此产生的误差较大，所以利用插值后的数据计算(10.13)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=l