Plot the integral of a discrete function

2 次查看(过去 30 天)
Good morning, I have to plot the integral of a discrete function (t,df/dt) defined in a file *.txt. I want to plot (t,f). f is an angle. I've tried with the command trapz but it gives only the numeric value of the area. Thanks

回答(2 个)

Youssef  Khmou
Youssef Khmou 2014-5-14
Trapz function returns a scalar value, numeric primitive can be calculated using the following function :
function itg=integral(f,dx,smooth);
% INTEGRAL ANS=INTEGRAL(F,DX)
% This function computes the integral of F(X)DX where the integrand
% is specified at discrete points F spaced DX apart (F is a vector,
% DX is a scalar). Simpsons Rule is used, so that the error
% is O(dx^5*F4). (F4 is the 4th derivative of F).
%
% If F is a matrix, then the integration is done for each column.
%
% If F is really spiky, then INTEGRAL(F,DX,'smooth') may
% provide a better looking result (the result is smoothed
% with a 3 point triangular filter).
%
% Author: RP (WHOI) 15/Aug/92
[N,M]=size(f);
if (N==1 | M==1),
N=max(size(f));
itg=zeros(size(f));
itg(1)=0; % first element
itg(2)=(5*f(1)+8*f(2)-f(3))*dx/12; % Parabolic approx to second
itg(3:N)=(f(1:N-2)+4*f(2:N-1)+f(3:N))*dx/3; % Simpsons rule for 2-segment
% intervals
itg(1:2:N)=cumsum(itg(1:2:N)); % Sum up 2-seg integrals
itg(2:2:N)=cumsum(itg(2:2:N));
if (nargin>2), % ... apply smoothing
itg(2:N-1)=(itg(1:N-2)+2*itg(2:N-1)+itg(3:N))/4;
itg(N)= (itg(N-1)+itg(N))/2;
end;
else
itg=zeros(size(f));
itg(1,:)=zeros(1,M);
itg(2,:)=(5*f(1,:)+8*f(2,:)-f(3,:))*dx/12;
itg(3:N,:)=(f(1:N-2,:)+4*f(2:N-1,:)+f(3:N,:))*dx/3;
itg(1:2:N,:)=cumsum(itg(1:2:N,:)); % Sum up 2-seg integrals
itg(2:2:N,:)=cumsum(itg(2:2:N,:));
if (nargin>2), % ... apply smoothing
itg(2:N-1,:)=(itg(1:N-2,:)+2*itg(2:N-1,:)+itg(3:N,:))/4;
itg(N,:)= (itg(N-1,:)+itg(N,:))/2;
end;
end;
Example :
t=0:0.1:10;
y=sin(t);
z=integral(y,0.1);
figure; plot(t,y,t,z);

Antonio
Antonio 2014-5-14
编辑:Antonio 2014-5-14
Thanks for the answer.
The code gives me the error
In an assignment A(I) = B, the number of elements in B and I must be the same.
could you apply your code to the real case of study I've attached here?
I need to plot f.
load data_hw1.txt
dx = data_hw1(:,1)
F = data_hw1(:,3) % F=df/dx
ps my homework suggests the trapezoidal integration.
  2 个评论
Youssef  Khmou
Youssef Khmou 2014-5-14
the sampling rate is :
DX=dx(300)-dx(299);
FF=integral(F,DX);
figure; plot(dx,F,dx,FF),legend(' F','\int F');
Antonio
Antonio 2014-5-14
编辑:Antonio 2014-5-14
It gives me error :( I've tried with this code (I need to plot angle in the interval 0deg:360deg. is it correct?
load data_hw1.txt
dx = data_hw1(:,1)
f = data_hw1(:,3);
FF=cumsum(f)
a=FF./10;
for n = 1:1:length(dx)
if a(n,1) < -360
FF(n)=(a(n,1)+720)
else
FF(n)=(a(n,1)+360)
end
end
figure(1)
plot(dx,f,dx,FF)

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by