Fourier Analysis on a Damped Harmonic Oscillator

16 views (last 30 days)
I need to test if my oscillator model is damped or not, and so I need to derive frequency in order to create a graph that provides me with the correct information. This is my code, I'm going to need to send all of it:
k = 1;
t = 0.01;
ind = 1;
tmax = 100;
v = 5;
gamma = 0.005;
h = 0.01;
t_array(1) = 0.01;
v_array(1) = 5;
x_pos_array(1) = 0;
while (t <= tmax)
a = -(2*gamma*v_array(ind))-(k*x_pos_array(ind));
x_pos_array(ind+1) =x_pos_array(ind) + v_array(ind)*h;
v_array(ind+1) = v_array(ind) + a*h;
t_array(ind+1) = t_array(ind) + h;
t=t_array(ind+1);
ind = ind + 1;
end
F = fft(x_pos_array);
AF = abs(F);
xF = x_pos_array;
for i = 1:10000
D(i) = AF(i)*t_array(i);
i = i + 1;
end
plot (D,x_pos_array)
Only x_pos_array affects the fourier transform, however the majority of it is so intertwined I had to send the whole thing.
Now, when comparing my results to known results (https://en.wikipedia.org/wiki/Harmonic_oscillator, https://au.mathworks.com/help/symbolic/physics-damped-harmonic-oscillator.html), the graph that this code spits out is nowhere similar to the graph that should be there.
Have I written the code relating to the fourier transformation and associated graph correctly? Or do I need to do something else?

Answers (1)

Mathieu NOE
Mathieu NOE on 14 Apr 2022
hello
I improved a bit the firts section of your code but I don't get what you are trying to do wit the second portion of the cde
F = fft(x_pos_array);
AF = abs(F);
xF = x_pos_array;
for i = 1:10000
D(i) = AF(i)*t_array(i);
i = i + 1;
end
plot (D,x_pos_array)
improved code :
clc
clearvars
dzeta = 0.05; % damping ratio
% physical constants
m = 10; % mass
k = 1; % stiffness
c = 2*dzeta*sqrt(k*m); % viscous damping
ind = 1;
tmax = 100;
v = 5;
h = 0.01;
t_array(1) = 0.01;
v_array(1) = 5;
x_pos_array(1) = 0;
t = t_array(1);
while (t <= tmax)
a = -1/m*((c*v_array(ind))+(k*x_pos_array(ind))); % acceleration = forces divided by mass
v_array(ind+1) = v_array(ind) + a*h; % update velocity first
x_pos_array(ind+1) = x_pos_array(ind) + v_array(ind+1)*h; % and after position
t_array(ind+1) = t_array(ind) + h;
t=t_array(ind+1);
ind = ind + 1;
end
plot(t_array,x_pos_array);

Categories

Find more on RF Toolbox in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!