# How to calculate average for a multiple of time period

8 views (last 30 days)
Abhik Saha on 24 Aug 2022
Commented: Mathieu NOE on 25 Aug 2022
Hello, Everyone, I have a second order diffrential equation from which I get position as a function of time and velocity as a function of time. (code below) But after that when I try to find the time period(T)=2*pi/omega from the data then I could not find a correct way to calculate that using the data generated. After that I want to calculate average of position and average of velocity. Now the average of position is define in this way
<x>=, same thing for the velocity where n1 is the integer number then I need to plot average value of x for different n1? I am new to matlab please help regarding this. I have copied the code below.
----------------main program----------------------------
t=linspace(0.2,10,1000);
y0=[1 0];
[tsol, ysol]=ode45(@(t,y0) firstodefun2(t,y0), t, y0);
figure(1);clf;
plot(tsol,ysol(:,1),'b')
ylabel('Position')
figure(2);clf;
plot(tsol,ysol(:,2),'r')
ylabel('velocity')
------------------function------------------------------
function dy=firstodefun2(t,y0)
G=1;gamma=1;omega=100;
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=G*sin(omega*t)-gamma*y0(2);
end

Mathieu NOE on 24 Aug 2022
here your are , my friend
the frequency is estimated by taking the fft of the velocity signal and getting the dominant peak.
but the average displacement definition is not correct : the integral must be divided by the time of integration to get the average (otherwise you get an integral of the displacement but not it's average value)
hope it helps
%----------------main program----------------------------
dt = 1e-3;
t=(0:dt:10);
samples = numel(t);
y0=[1 0];
[tsol, ysol]=ode45(@(t,y0) firstodefun2(t,y0), t, y0);
pos = ysol(:,1);
velo = ysol(:,2);
figure(1);clf;
plot(tsol,pos,'b')
ylabel('Position')
figure(2);clf;
plot(tsol,velo,'r')
ylabel('velocity')
% get omega back from velocity fft
Fs = 1/dt;
[f,P] = one_sidded_fft(velo,Fs);
% find the peak frequency
[PKS,LOCS] = findpeaks(P,'MinPeakHeight',max(P)/2);
freq = f(LOCS);
omega_fft = 2*pi*freq
figure(3),
plot(f,P)
title('velocity fft')
xlabel('Frequency (Hz)');
ylabel('Amplitude')
% average Position (integral over n1 samples from t0)
n1 = 100; % must be less or equal to samples
t1 = t(1:n1);
pos1 = pos(1:n1);
av_Position = trapz(t1,pos1)/t1(end)
%------------------function------------------------------
function dy=firstodefun2(t,y0)
G=1;gamma=1;omega=100;
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=G*sin(omega*t)-gamma*y0(2);
end
function [f,P] = one_sidded_fft(X,Fs)
samples = length(X);
% regular (non windowed) fft%
P = abs(fft(X)/samples);
% %% hanning window fft%
% P = abs(fft(X.*hanning(samples))*4/samples);
% one sidded fft spectrum % Select first half
if rem(samples,2) % nfft odd
select = (1:(samples+1)/2)';
else
select = (1:samples/2+1)';
end
P = P(select);
P(2:end-1) = 2*P(2:end-1);
f = (select - 1)*Fs/samples;
end
##### 2 CommentsShow NoneHide None
Abhik Saha on 25 Aug 2022
Thank you for the code and the justification. yes you are right to calculate average we need to divide by time.
Thanks!
Mathieu NOE on 25 Aug 2022
My pleasure !

### Categories

Find more on Mathematics in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!