필터 지우기
필터 지우기

Plot the integral of a discrete function

조회 수: 2 (최근 30일)
Antonio
Antonio 2014년 5월 14일
편집: Antonio 2014년 5월 14일
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