Un metodo: RE:ugente integrar
Hola Joakin
Te envio un codigo que encontre en la web, podrias probar con ell es un metodo de integracion por Simpson o revizar en la pagina web de Matlab:
http://www.mathworks.com/matlabcentral/fileexchange/loadCategory.do?objectType=category&objectId=16&objectName=Integration
Saludos
Fismat
--------------------------------
function z = simps(x,y,dim)
%SIMPS Simpson's numerical integration.
% Z = SIMPS(Y) computes an approximation of the integral of Y via the
% Simpson's method (with unit spacing). To compute the integral for
% spacing different from one, multiply Z by the spacing increment.
%
% For vectors, SIMPS(Y) is the integral of Y. For matrices, SIMPS(Y) is a
% row vector with the integral over each column. For N-D arrays, SIMPS(Y)
% works across the first non-singleton dimension.
%
% Z = SIMPS(X,Y) computes the integral of Y with respect to X using the
% Simpson's rule. X and Y must be vectors of the same length, or X must
% be a column vector and Y an array whose first non-singleton dimension
% is length(X). SIMPS operates along this dimension.
%
% Z = SIMPS(X,Y,DIM) or SIMPS(Y,DIM) integrates across dimension DIM of
% Y. The length of X must be the same as size(Y,DIM).
%
% Example:
% -------
% If X = [0 1 3 6] then simps(X,X.^2) is 72 whereas trapz(X,X.^2) is 78
%
% If Y = [0 1 2
% 3 4 5
% 6 7 8]
% then simps(Y,1) is [6 8 10] and simps(Y,2) is [2; 8; 14]
%
% Class support for inputs X, Y:
% float: double, single
%
% Damien Garcia, 08/2007, directly adapted from TRAPZ
%
% See also CUMSIMPS, TRAPZ, QUAD.
%% Make sure x and y are column vectors, or y is a matrix.
perm = []; nshifts = 0;
if nargin == 3 % simps(x,y,dim)
perm = [dim:max(ndims(y),dim) 1:dim-1];
y = permute(y,perm);
[m,n] = size(y);
elseif nargin==2 && isscalar(y) % simps(y,dim)
dim = y; y = x;
perm = [dim:max(ndims(y),dim) 1:dim-1];
y = permute(y,perm);
[m,n] = size(y);
x = 1:m;
else % simps(y) or simps(x,y)
if nargin < 2, y = x; end
[y,nshifts] = shiftdim(y);
[m,n] = size(y);
if nargin < 2, x = 1:m; end
end
x = x(:);
if length(x) ~= m
if isempty(perm) % dim argument not given
error('MATLAB:simps:LengthXmismatchY',...
'LENGTH(X) must equal the length of the first non-singleton dimension of Y.');
else
error('MATLAB:simps:LengthXmismatchY',...
'LENGTH(X) must equal the length of the DIM''th dimension of Y.');
end
end
if m<3
error('MATLAB:simps:LengthX',...
'LENGTH(X) must be at least 3.');
end
%% The output size for [] is a special case when DIM is not given.
if isempty(perm) && isequal(y,[])
z = zeros(1,class(y));
return;
end
%% Simpson's rule
dx = repmat(diff(x,1,1),1,n);
dx1 = dx(1:end-1,:);
dx2 = dx(2:end,:);
alpha = (dx1+dx2)./dx1/6;
a0 = alpha.*(2*dx1-dx2);
a1 = alpha.*(dx1+dx2).^2./dx2;
a2 = alpha.*dx1./dx2.*(2*dx2-dx1);
z = sum(a0(1:2:end,:).*y(1:2:m-2,:) +...
a1(1:2:end,:).*y(2:2:m-1,:) +...
a2(1:2:end,:).*y(3:2:m,:),1);
if rem(m,2) == 0 % Adjusting if length(x) is even
C = vander(x(end-2:end))\y(end-2:end,:);
z = z + C(1,:).*(x(end,:).^3-x(end-1,:).^3)/3 +...
C(2,:).*(x(end,:).^2-x(end-1,:).^2)/2 +...
C(3,:).*dx(end,:);
end
%% Resizing
siz = size(y); siz(1) = 1;
z = reshape(z,[ones(1,nshifts),siz]);
if ~isempty(perm), z = ipermute(z,perm); end