《数值积分的matlab实现》由会员分享，可在线阅读，更多相关《数值积分的matlab实现（13页珍藏版）》请在金锄头文库上搜索。
1、实验10数值积分实验目的:1了解数值积分的基本原理；2 .熟练掌握数值积分的 MATLAB实现；3 会用数值积分方法解决一些实际问题。实验内容：积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。同微分一样，在微积分中，它也是通过极限定义的， 由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所1 sin x以仍然得不到积分的精确值，如不定积分dx。这时我们一般考虑用数值方法计算其0 x近似值，称为数值积分。10.1 数值微分简介设函数y f (x)在x*可导，则其导数为XXfmHhhX(10.1如果函数y f(
2、x)以列表形式给出(见表10-1 ),则其精确值无法求得，但可由下式求得其近似值*f (x h) f (x )f(X )(10.2 )h表 10-1xx0x1X?xny0y1y2Vny般的，步长h越小，所得结果越精确。(10.2 )式右端项的分子称为函数y f (x)在x的差分，分母称为自变量在x的差分，所以右端项又称为 差商。数值微分即用差商近似 代替微商。常用的差商公式为:f (Xo)f(x h) f (Xo h)2hf (Xo)3y 4yi y22hf (Xn)yn 2 4yn 1 3yn2h(10.3 )(10.4 )(10.5 )其误差均为0(h )，称为统称三点公式。10.2 数值
3、微分的MATLAB实现MATLAB提供了一个指令求解一阶向前差分，其使用格式为:dx=diff(x)其中x是n维数组，dx为n 1维数组x2为，怡X2,，Xn为，这样基于两点的数值导数可通过指令diff(x)/h 实现。对于三点公式，读者可参考例 1的M函数文件diff3.m。例1用三点公式计算 y f (x)在x 1.0,121.4处的导数值，f(x)的值由下表给出。Xf(x)解：建立三点公式的M函数文件diff3.m 如下:fun cti on f=diff3(x,y)n=le ngth(x);h=x(2)-x(1); f(1)=(-3*y(1)+4*y (2)-y(3)/(2*h);fo
4、r j=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指令窗中输入指令:运行得各点的导数值为： ，。所以 y f (x) 在 x ，和。对于高阶导数， MATLAB 提供了几个指令借助于样条函数进行求导，详细使用步骤如 下：step1 ：对给定数据点( x,y )，利用指令 pp=spline(x,y) ，获得三次样条函数数据 pp ， 供后面 ppval 等指令使用。其中， pp 是一个分段多项式所对应的行向量，它包含此多项式 的阶数、段数、节点的横坐标值和各段多项式
5、的系数。step2 ：对于上面所求的数据向量 pp ，利用指令 breaks,coefs,m,n=unmkpp(pp) 进 行处理，生成几个有序的分段多项式 pp 。step3 ：对各个分段多项式 pp 的系数，利用函数 ppval 生成其相应导数分段多项式的 系数，再利用指令 mkpp 生成相应的导数分段多项式step4 ：将待求点 xx 代入此导数多项式，即得样条导数值。上述过程可建立 M 函数文件 ppd.m 实现如下：function dy=ppd(pp)breaks,coefs,m=unmkpp(pp);for i=1:mcoefsm(i,:)=polyder(coefs(i,:);
6、enddy=mkpp(breaks,coefsm);于是，如果已知节点处的值 x,y ，可用下面指令计算 xx 处的导数 dyy ：pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx);例 2 基于正弦函数 y sin x 的数据点，利用三点公式和三次样条插值分别求导，并与解析所求得的导数进行比较。解:编写M脚本文件bijiao.m 如下：h=0.1*pi;x=0:h:2*pi;y=si n(x);dy仁diff3(x,y);pp=spli ne(x,y);dy=ppd(pp);dy2=ppval(dy,x);z=cos(x);error1= no rm(dy1
7、-z),error2=n orm(dy2-z)plot(x,dy1,k:,x,dy2,r-,x,z,b)运行得结果为：，生成图形见图10.1 o图10.1三点公式、三次样条插值与解析求导比较图显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图2.15可知利用三次样条插值求导所得曲线与解析求导曲线基本重合，而三点公式在极值点附近和两个端点附近误差较大，其它点吻合的较好。10.3 应用示例：湖水温度变化问题问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的 温度越低。这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。 对某个湖的水温进行观测得
8、数据见表10-2 o表10-2某湖的水温观测数据深度(m )0温度C)试找出湖水温度变化最大的深度。1 问题的分析湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导对于给定的数据,数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。可以利用数值微分计算各深度的温度变化值， 从而得到温度变化最大的深度， 但考虑到所给 的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值， 以得到更准确的结果。2模型的建立及求解记湖水的深度为h (m )，相应的温度为T (C),且有T T(h)，并假定函数T(h)可 导。对给定的数据进行三次样条插值
9、，并对其求导，得到T(h)的插值导函数；然后将给定的深度数据加密， 搜索加密数据的导数值的绝对值， 找出其最大值及其相应的深度， 相应的 MATLAB 指令如下：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=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.),grid on运行得导函数绝对值的最大值点为： h
10、=11.4 ，最大值为 1.6139,即湖水在深度为时温 度变化最大，如图 10.2 所示(黑点为温度变化最大的点) 。图 10.2 湖水温度变化曲线图10.4 数值积分简介考虑定积分bf (x)dx ( 10.6 )a如果被积函数f(x)是以列表形式给出，则其求解思想同数值微分类似，即用逼近多项b式Pn(x)近似地代替被积函数f (x)，然后计算积分R(x)dx，得(10.6 )式的近似值；a如果被积函数的原函数不是初等函数， 则将积分区间进行细分， 对每个小区间， 用一个近似函数代替被积函数 f(x)，然后积分得(10.6 )式的近似值。这两种类型最终都可归结为函 数f (x)在节点 兀上
11、的函数值f(xj的某种线性组合，即下面数值求积公式：f(x)dxaAkf(xk)k 0(107)f(x)dxAJ (xQk 0(10.8)其中R f为截断误差。此误差可用代数精度衡量，代数精度越高，误差越小；反之误差越大。代数精度是用来衡量数值积分公式近似程度的办法，如果f (x)是一个次数不超过 m的代数多项式，(10.7 )式等号成立；而当 f (x)是一个m 1次多项式时，(10.7 )式不能精 确成立，则称(10.7 )式的代数精度为 m。选取不同的近似函数，可产生不同的数值求积公式，常见的有：梯形公式、辛普森公式 和高斯公式。10.5 数值积分的MATLAB 实现MATLAB提供了下
12、面几个函数计算积分，其使用格式分别为：(1) trapz(x)采用梯形公式计算积分(h 1),x为fk(k 0,1,n)(2) quad(fun,a,b,tol)采用自适应 Simpson 法计算积分(3) quadl(fun,a,b,tol)采用自适应 Gauss-Lobatto 法计算积分其中fun为被积函数；tol是可选项，表示绝对误差，a，b为积分的上、下限。1 2例1分别利用梯形公式、Simpson 公式和Gauss-Lobatto 法计算 0 + 1 x dx，并与 其精确值比较。解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方法求 得的数值解比较，其 MAT
13、LAB指令为:syms xxzO=simple(i nt(sqrt(1+xxA2),0,1)z=double(z0);z=vpa(z,8)x=0:0.01:1;y=sqrt(1+x.A2);z1=trapz(y)*0.01;z 1=vpa(z1,8),err1=z-z1;err1=vpa(err1,8)z2=quad(sqrt(1+x.A2), 0,1);z2=vpa(z2,8),err2=z-z2;err2=vpa(err2,8)z3=quadl(sqrt(1+x.A2), 0,1);z3=vpa(z3,8),err3=z-z3;err3=vpa(err3,8)运行得精确值为*C，2In(一
14、21)，三种公式计算得数值积分值分别为，和，其相应误差分别为，.1e-6和0.，由三者误差可见，Gauss-Lobatto法计算最为精确，Simpson 公式次之，梯形公式最差，但它也能精确到小数点后5位数。例2人造地球卫星轨道可视为平面上的椭圆。我国第一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为 6371km，求该卫星的轨道长度。解：卫星轨道椭圆的参数方程为 x a cost, y bsin t(0 t 2 ), a,b分别是长、短半轴，则根据所给数据知a =6371+2384=8755, b =6371+439=6810。由对弧长的曲线积分知识知，椭圆的长度为1L 4 02 (a2 sin2t b2cos2 t)2dt上积分称为椭圆积分，它无法用解析方法计算，可用计算其数值解，编写M函数文件如下：fun cti on y=y(t)a=8755;b=6810;y=4*sqrt(aA2*si n( t).A2+bA2*cos(t).A2);在MATLAB指令窗中输入以下指令：l=quad(y,0,pi/2)对于用列表形式给出的函数，上述方法不再适用，可利用指令spli ne构造三次样条插值函数，再计算积分，具体步骤可参考例2。例3在桥梁的一端每隔一段时间记录1min有几辆车过桥，得到数据见表10-3 :表10-3过桥车辆数据
