Fourier Analysis on a Damped Harmonic Oscillator

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?

Answer 1

[Mathieu NOE]

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);