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

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 rowx=x( : ); nx=length(x); % make x a columnry , cy=size(y);if ry=1&cy=nx , y=y . ; ry=cy ; cy=1 ; flag=1 ; endif nx=ry , error( X and Y not the right size ) , enddx=x(2 : nx)-x(1 : nx-1); % width of each trapezoiddx=dx( : , ones(1 , cy); % duplicate for each column in yyave=(y(2 : ry , : )+y(1 : ry-1 , :)/2; % average of heightsz=zeros(1 , cy); cumsum(dx .* yave); % Use cumsum to find areaif 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 the derivative 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.Weighted central difference are used.% Quadratic approximation is used at the endpoints.% % See also mmintgrl% Copyrigth (c) 1996 by Prentice-Hall , Inc.flag=0;% flag is True if y is a rowx=x( : ); nx=length(x); % make x a columnry , cy=size(y);if ry=1&cy=nx , y=y .; ry=cy; cy=1; flag=1; endif nx=ry , error( X and Y not the right size ) , endif nx=z1;for i=1 : cy % fit quadratic at endpoints of each column p1=polyfit(x(1 : 3) , y(1 : 3 , i) , 2); % quadratic at first point z1(i)=2*p1(1)*x(1)+p1(2); % evalute poly derivative pn=polyfit(x(nx-2 : nx) , y(ry-2 : ry , i) , 2); % quadratic at last point zn(i)=2*pn(1)*x(nx)+pn(2); % evaluate poly derivativeendz=z1; z; zn;if flag , z=z; end % if y was a row , return a row最后，给出一个例子：x=linspace(0 , 2*pi , 30)y=sin(x); % create datayi=mmintgrl(x , y); % find integralyd=mmderiv(x , y); % find derivativeplot(x , y , x , yi , - , x , yd , : ) % plot results注意这个积分定性地证明了等式：而微分定性地证明了等式：图13.11 y=sin(x)极其积分、微分曲线13.8 小结表13.1总结了本章所讨论的函数。表13.1数值分析函数fplot( fname , lb ub)绘出上下限之间的函数fmin( fname , lb ub)寻找上下限内的标量最小值fimis( fname , xo)寻找xo附近的向量最小值fzero( fname , xo)寻找x