trapc.m

来自「用于求积分」· M 代码 · 共 17 行

function [Itcp]=trapc(a,b,m,f,choice,varargin)
h=(b-a)/m;
switch choice
    case 1,
    x=linspace(a,b,m+1);
    fpm=feval(f,x,varargin{:});
    fpm(2:end-1)=2*fpm(2:end-1);
    Itpc=0.5*h*sum(fpm);
    case 2,
        z=linspace(a,b,m+1);
        x=z(1:end-1)+h/2*(1-1/sqrt(3));
        x=[x,x+h/sqrt(3)];
        fpm=feval(f,x,varargin{:});
        Itpc=0.5*h*sum(fpm);
    otherwise
        disp('unknown formula');
end