積分用梯形规则计算。用加权中心差分计算微分及数值分析函数.docVIP

13.7? M文件举例 这里所介绍的《精通MATLAB工具箱》中的M文件可近似求解由采样值给出的函数的积分和微分。这里假定这些函数本身不存在，且独立变量也许不是线性间隔。例如，已装载到MATLAB中要分析的数据来源于实验测试。 对于所包含的数据缺乏函数描述，有许多种积分和微分的方法。如前所述，人们可以用最小二乘多项式拟合数据，然后在多项式的描述上进行操作。另一种方法是寻找数据的三次样条表示，然后运用《精通MATLAB工具箱》中的函数spintgrl和spderiv来分别寻找积分和微分的样条表示。这里所介绍的方法提供了另一种更简单的方法。积分用梯形规则计算。用加权中心差分计算微分。此外，将函数设计成在矩阵形式下工作，矩阵的列代表各与自变量有关的因变量。 正如这章前面所述，MATLAB函数trapz计算在某有限区间的梯形积分。这里我们寻找的积分是自变量为x的函数。即如果y=f(x)，我们寻找： 式中的x1是向量x的第一个元素。用梯形规则，这个积分近似为： ? 且S(x1)=0 这样，第k个数据点的积分是上述梯形面积的累加和。函数mmintgrl实现的这个算法如下： function z=mmintgrl(x , y) %? MMINTgrl Compute Integral using Trapezoidal Rule. %?? MMINTGRL(X , Y) computes the integral of the function y=f(x) given the %?? data in X and Y. X must be a vector ,? but Y may be a column oriented %?? data matrix. The length of X must equal the length of Y if Y is a %?? vector ,? or it must equal the number of rows in Y if Y is a matrix. %? %?? X need not be equally spaced. The trapezoidal algorithm is used. %? %?? See also mmderiv %?? Copyrigth (c) 1996 by Prentice-Hall , Inc. flag=0; %?? falg is True if y is a row x=x( : ); nx=length(x);? %?? make x a column [ry , cy]=size(y); if ry==1cy==nx , y=y . ; ry=cy ; cy=1 ; flag=1 ; end if nx~=ry ,? error( X and Y not the right size ) , end dx=x(2 : nx)-x(1 : nx-1);???????????? %?? width of each trapezoid dx=dx( : , ones(1 , cy));????????????? %?? duplicate for each column in y yave=(y(2 : ry , : )+y(1 : ry-1 , :))/2;???????????????? %?? average of heights z=[zeros(1 , cy); cumsum(dx .* yave)];??????????? %?? Use cumsum to find area if flag , z=z; end %?? if y was a row , return a row 在介绍上述函数的使用之前，考虑微分。在这种情况下，人们感兴趣的就是刚给定数据点的近似斜率。这里介绍一种下述的中心差分 图13.10 加权中心差分方法 从图13.11可知，在第k个点的近似微分是： ?????? 式中 ， 并且Mk是连接yk-1到yk的直线的斜率。这样，第k点的微分是相邻两点间斜率的加权平均，离该点越近的点权越重。在第一个和最后一个数据点上，不能简单按照上述方法进行处理，因为这两个点都没有伴随的直线段。对于这些数据点，需要用另外的方法。这里所采取的方法是用二次多项式拟合前3个点（或最后3个点），并且计算这个多项式第一个（或最后一个）点的微分。函数mmderiv实现的这个算法如下： function z=mmderiv(x , y) %? MMDERIV Compute Derivative Using Weighted Central Differences. %?? MMDERIV(X , Y) computes