Plot the integral of a discrete function

3 views (last 30 days)
Antonio
Antonio on 14 May 2014
Edited: Antonio on 14 May 2014
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

Answers (2)

Youssef  Khmou
Youssef Khmou on 14 May 2014
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 on 14 May 2014
Edited: Antonio on 14 May 2014
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 Comments
Youssef  Khmou
Youssef Khmou on 14 May 2014
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 on 14 May 2014
Edited: Antonio on 14 May 2014
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)

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!