% Square Wave problem
% INPUTS
T = 5; % sec, Period of forcing function
F_o = 30; % N, amplitude
k = 1000; % N/m, spring stiffness
omega_n = 3; %rad/s, natural frequency
zeta = 0.1; % damping factor
num_terms = 5; % Number of terms for Fourier series
num_periods = 2;
del_t = T/100; % Time increment for results
%% CALCULATIONS
time = 0:del_t:num_periods*T;
omega_T = (2*pi) / T;
% Exact Forcing Function
for i = 1:length(time)
if time(i) <= T/2
F(i) = -F_o;
elseif time(i) <= T
F(i) = F_o;
elseif time(i) <= (3*T)/2
F(i) = -F_o;
else
F(i) = F_o;
end
end
% Fourier approximation of Forcing Function
for n = 1:num_terms
aj = zeros(n,1); % numerical integration
bj = [-120/pi; 0; -40/pi; 0; -24/pi];
end
for i = 1:length(time)
F_1(i) = aj(1)*cos(1*omega_T*time(i)) + bj(1)*sin(1*omega_T*time(i)) ;
F_2(i) = aj(2)*cos(2*omega_T*time(i)) + bj(2)*sin(2*omega_T*time(i)) ;
F_3(i) = aj(3)*cos(3*omega_T*time(i)) + bj(3)*sin(3*omega_T*time(i)) ;
F_4(i) = aj(4)*cos(4*omega_T*time(i)) + bj(4)*sin(4*omega_T*time(i)) ;
F_5(i) = aj(5)*cos(5*omega_T*time(i)) + bj(5)*sin(5*omega_T*time(i)) ;
F_Fourier = F_1 + F_2 + F_3 + F_4 + F_5;
end
% Response
r = zeros(1,5);
tan = zeros(1,5);
alpha = zeros(1,5);
x_1 = zeros(1,201);
x_3 = zeros(1,201);
x_5 = zeros(1,201);
for s = 1:num_terms
r(s) = s.*(omega_T / omega_n);
tan(s) = (2*zeta.*r) / (1-(r.^2));
alpha(s) = atan(tan(s));
end
X_1 = (bj(1) / k) * (1 / sqrt( ((1 - (r(1)^2))^2) + ((2*zeta*(r(1)))^2)));
% X_2 = (bj(2) / k) * (1 / sqrt( ((1 - (r(2)^2))^2) + ((2*zeta*(r(2)))^2)));
X_3 = (bj(3) / k) * (1 / sqrt( ((1 - (r(3)^2))^2) + ((2*zeta*(r(3)))^2)));
% X_4 = (bj(4) / k) * (1 / sqrt( ((1 - (r(4)^2))^2) + ((2*zeta*(r(4)))^2)));
X_5 = (bj(5) / k) * (1 / sqrt( ((1 - (r(5)^2))^2) + ((2*zeta*(r(5)))^2)));
x_1 = X_1*(sin(((s(1)*omega_T*time) - alpha(s(1)))));
% x_2 = X_2*(sin(((s(2)*omega_T*time) - alpha(s(2)))));
x_3 = X_3*sin(((s(3)*omega_T*time) - alpha(s(3))));
% x_4 = X_4*sin(((s(4)*omega_T*time) - alpha(s(4))));
x_5 = X_5*sin(((s(5)*omega_T*time) - alpha(s(5))));
x_p = x_1 + x_3 + x_5;
