Question 4.7: Represent the triangular function shown in Figure 4.12 in te...
Represent the triangular function shown in Figure 4.12 in terms of Fourier series using MATLAB.
The MATLAB code is written in the following parts: Main module DOS4_1 calls MATLAB function ‘quad’. The function ‘quad’ numerically evaluates an integral using Simpson rule.
Function curve1 calculates the constant term a_{0}
Function curve1c calculates the a_{n} coefficients of cosine terms
Function curve1s calculates the b_{n} coefficients of sine terms
The Fourier representation with 15 terms is shown in Figure 4.13.
------------------
CODE DOS4_1 Triangular Function
------------------
% Dynamics of Structures by Ashok K. Jain I. I. T. Roorkee
% DOS4_1
% Fourier representation of a periodic force: Triangular function
% Calls: curve1, curve1c and curve1s
% curve1c: cosine constants; curve1s: sine constants
% The program calls function “quad” to integrate and find Fourier constants
% N: total number of terms to be included
% T: period of periodic force
% P: Amplitude of periodic force
%
close all;
clear all;
clc;
P = 10;
T = 2;
fprintf(‘Amplitude of periodic force = %.3g N\n’,P);
fprintf(‘Period of force = %.3g sec\n’,T);
a0 = quad(@(t)curve1(t,T,P),0,T,1e-9);
% for N=15:15
for N=1:10
for n=1:N
a(n) = quad(@(t)curve1c(t,n,T,P),0,T,1e-9);
b(n) = quad(@(t)curve1s(t,n,T,P),0,T,1e-9);
end
Outp = 0;
x=1;
for t=0:0.05:T*2
c=0;
s=0;
for n=1:N
c=c+a(n)*cos(t*2*pi*n/T);
s=s+b(n)*sin(t*2*pi*n/T);
end
Outp(x)=s+c+a0;
x=x+1;
end
subplot (5,2,N)
plot(0:0.05:T*2,Outp,’.-’)
grid on
pause(1);
end
------------------
CODE CURVE1
------------------
% curve1
% Called by: DOS4_1
% Fourier representation of a periodic force: Triangular function
%
function a = curve1(t,T,P)
w = 2*pi/T;
for d=1:length(t)
if t(d)<=T/2
p1(d) = 2*t(d)*P/T;
else
p1(d) = 2*(T-t(d))*P/T;
end
end
a = (1/T)*p1;
end
------------------
CODE CURVE1C
------------------
% curve1c
% Called by: DOS4_1
% Fourier representation of a periodic force: Triangular function
% Computes constants for cosine terms
%
function a = curve1c(t,n,T,P)
w = 2*pi/T;
for d=1:length(t)
if t(d)<=T/2
p1(d) = 2*t(d)*P/T;
else
p1(d) = 2*(T-t(d))*P/T;
end
end
a = (2/T)*p1.*cos(w*t*n);
end
------------------
CODE CURVE1S
------------------
% curve1s
% Called by: DOS4_1
% Fourier representation of a periodic force: triangular function
% Computes constants for sine terms
%
function b = curve1s(t,n,T,P)
w = 2*pi/T;
for d=1:length(t)
if t(d)<=T/2
p1(d) = 2*t(d)*P/T;
else
p1(d) = 2*(T-t(d))*P/T;
end
end
b = (2/T)*p1.*sin(w*t*n);
end