数值积分的matlab实现

实验10数值积分实验目的：1 .理解数值积分的基本原理2 .熟练运用数值积分的MATLAB的实现3 .用数值积分法解决一些实际问题。实验内容：积分是数学的基本概念，也广泛应用于实际问题。 与导数类似，在微积分也是限定性的，然而不能直接对实际问题中遇到的函数进行积分，因为通常是以列表形式给出的。 另外，虽然有解析式，但是因为原函数不是初等函数，所以也有不能得到不定积分这样的积分的正确值的函数。 在这种情况下，一般考虑到以数值方式计算近似值，将其称为数值积分。10.1数值微分概要如果函数是导电的，则该导数是(10.1 )函数以列表形式给出时(参照表10-1 )，无法求出其正确的值，但其近似值可以通过下式求出(10.2 )表10-1一般来说，步长越小，结果越准确。 (10.2 )式的右端项的分子称为函数的差分，分母称为参数的差分，因此将右端项称为差分商。 用差商代替微分商近似数值微分。 一般差饷公式如下：(10.3 )(10.4 )(10.5 )其误差全部统称为三点式。10.2数值微分的MATLAB实现MATLAB提供了使用以下格式确定主前向差异的指令dx=diff(x )这里，x是三维数组，dx是三维数组，基于两个点的数值导数可以由命令diff(x)/h来实现。 对于三点式，读取者可参考例1的m函数文件diff3.m。例13点式，计算1.0、1.2、1.4下的微分值，其值如下表所示。1.01.1.2.3.40.2500 0.2268 0.2066 0.1890 0.1736解：作成三点式的m函数文件diff3.m如下所示函数f=diff3 (x，y )n=长度(x ) h=x (2)-x (1)f(1)=(-3*y(1) 4*y(2)-y(3)/(2*h )for j=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 )在“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。 因此，在1.0、1.2、1.4下导数值分别为-0.2470、-0.1890、-0.0014 .在MATLAB中，可以使用样条函数为高阶导数提供一些命令。 详细步骤如下：step1:对于数据点(x，y )，使用指令pp=spline(x，y )，得到三次样条函数数据pp，用于以后的ppval等指令。 在此，pp是与段多项式对应的行向量，包含该多项式的次数、段数、节点的横轴值、各段多项式的系数。步骤2 :对于以上求出的数据向量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:m米coefsm(i，)=polyder(coefs(i， ) )结束dy=mkpp(breaks，coefsm )并且，如果知道节点的值x、y，则可以用以下的指令计算xx的导数dyypp=spline(x，y )，dy=ppd(pp) dyy=ppval(dy，xx )例2 :通过三点式和三次样条插值分别获得基于正弦函数的数据点，并将其与通过分析获得的导数进行比较。在：中创建的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.718.322.927.2温度()22.822.822.820.613.911.711.111.1让我们来寻找湖水温度变化最大的深度。1 .问题分析由于湖水温度可视为与深度有关的函数，湖水温度的变化问题成为温度函数的导数问题，证明与导数的最大绝对值对应的深度是温度变化的最大深度。 对于给定数据，可以利用数值导数来计算每个深度的温度变化值，但是，考虑到给定数据的量小，通过插值方式来计算加密的深度数据的导数值，以获得更准确的结果，因为所计算的深度不准确。2 .建立和解决模型假定湖水的深度是(m )，相应的温度是()，并且函数得到。对给定数据执行三次样条插值，确定三次样条插值，对所获得的插值导数以及给定深度数据进行加密，搜索经加密数据的导数的绝对值，并找到其最大值和相应深度。 对应的MATLAB指令如下h=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2； t= 22.822.822.820.613.911.7.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.)，grid on取导数绝对值的最大值点为=11.4，最大值为1.6139，即湖水在深度11.4m处温度变化最大，如图10.2所示(黑点为温度变化最大的点)。图10.2湖水温度变化图10.4数值积分概要考虑定积分(10.6 )如果以列表形式给出了被积分函数，则其思想与数值微分类似，即将近似多项式近似地替换为被积分函数来计算积分，如果得到(10.6 )式的近似值的积分函数的原函数不是初等函数，则将积分区间细分，在各单元之间将积分函数替换为近似函数，从而得到(10.6 )式的近似值两种类型最终都可以导致函数的节点上的函数值的线性组合，即下一个数值乘积公式或(10.7 )(10.8 )这里是截止误差。 这个误差可以用代数精度来测量，代数精度越高误差越小，相反误差越大。代数精度是测量数值积分式近似程度的方法，如果是不超过次数的代数多项式，则(10.7 )式的等号成立，另一方面，如果是次多项式，则(10.7 )式不正确成立时，称为(10.7 )式的代数精度。如果选择不同的近似函数，则会生成不同的数值乘积表达式，这些表达式常见于梯形、辛普森表达式和高斯表达式。10.5数值积分的MATLAB实现MATLAB具有使用以下函数计算积分的功能(1)trapz(x )以梯形式计算积分()，x为(2)quad(fun，a，b，tol )使用适应辛普森法计算积分(3)使用3)quadl(fun，a，b，tol )自适应Gauss-Lobatto法计算积分fun是可选的乘积函数tol，表示绝对误差，a、b是积分的上下限。例1分别用梯形式、辛普森式、Gauss-Lobatto法计算，并与其正确的值进行了比较。解：对积分进行编码运算，将其计算结果转换为数值型，然后与用这3种方法求出的数值解进行比较，MATLAB命令如下syms xxz0=simple(int(sqrt(1 xx2)，0，1 ) )z=双精度(z0 ) z=VPA (z，8 )x=0:0.01:1； y=sqrt(1 x.2)z1=trapz(y)*0.01； Z1=VPA (Z1，8，8 )，err1=z-z1； err1=VPA (err 1，8 )z2=quad(sqrt(1 x.2)，0，1 ) z2=VPA (z2，8 )，err2=z-z2； err2=VPA (err 2，8 )z3=quadl(sqrt(1 x.2)，0，1 ) z3=VPA (z3，8 )，err3=z-z3； err3=VPA (err 3，8 )正确的运行值为1.1477936，用3个公式计算的数值积分值分别为1.1477995、1.1477935和1.1477936，其对应误差分别为-.59e-5、1e-6和0 .从3个误差可以看出，用Gauss-Lobatto法计算是最正确的例2人造地球卫星轨道可视为平面上的椭圆。 我国第一颗人造地球卫星的近地点距地球表面439km，远地点距地球表面2384km，地球半径6371km，求出了该卫星的轨道长度。解：卫星轨道椭圆的参数方程分别为长轴和短轴，并且根据给定数据=6371 2384=8755，=6371 439=6810。根据对弧长曲线积分的知识，椭圆长度上积分称为椭圆积分，无法用分析方法计算，因此计算数值解，并按如下方式创建m函数文件函数y=y (t )a=8755； b=6810；y=4* sqrt (a 2* sin (t ) 2b 2* cos (t ) 2；在MATLAB命令窗口中，输入以下命令l=quad(y，0，pi/2 )执行结果是l=4.9090e 004。 即卫星轨道长度为49090km。对于以列表形式提供的函数，上述方法将不应用，因此可以利用指令spline构建三次样条插值函数以计算积分，将参考实例2描述特定过程。例3记录桥梁的一端每小时有几辆车过桥，得到数据的是表10-3表10-3过桥车辆数据时间0:002:004:005:006:007:008:00车辆数/车辆22025825时间9:0010:3011:3012:3014:0016:0017:00车辆数/车辆12510127928时间18:0019:0020:0021:0022:0023:0024:00车辆数/车辆2210911893估计一天通过桥的车流量。解：记录时刻为时，如果该车辆数为，则求出的车流量为累计点，在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为：函数VF=函数(x，pp )vf=ppv