Fourier transform of symbolic function

13 views (last 30 days)
Hi, I want to conduct a frequenct analysis using the following code:
clear
syms t
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
k = 4.5682;
m = 0.2;
w = 0.0014;
J = 0.0327;
I = 0.00362;
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
fplot(fft(qwe.y))
fplot(fft(qwe.z))
However, when I am using fft(), I need to use numerical data. Is there a way to conduct a symbolic fourier transform, or do I need to convert the symbolic function to numeric data?
If I do need to convert my symbolic function to numeric data, how do I do this?

Accepted Answer

Star Strider
Star Strider on 3 Nov 2021
Use the fourier function. I’m guessing as to what ‘wilberforce_y’ and ‘wilberforce_z’ are.
syms t omega T
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
k = 4.5682;
m = 0.2;
w = 0.0014;
J = 0.0327;
I = 0.00362;
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
wilberforce_y = simplify(qwe.y, 500)
wilberforce_y = 
Fy = simplify(fourier(wilberforce_y), 500)
Fy = 
wilberforce_z = simplify(qwe.z, 500)
wilberforce_z = 
Fz = simplify(fourier(wilberforce_z), 500)
Fz = 
figure
fplot(wilberforce_y, [0 25])
grid
figure
fplot(fourier(wilberforce_y), [-10*pi 10*pi])
figure
fplot(wilberforce_z, [0 25])
grid
figure
fplot(fourier(wilberforce_z), [-10*pi, 10*pi])
That they are plotted as straight lines at 0 are likely due to the δ funcitons.
Integrating them manually as:
F2y = int(wilberforce_y * exp(1j*omega*t), t, -T, T)
F2y = 
F2y = vpa(simplify(F2y, 500), 5)
F2y = 
F2y = subs(F2y, T, 1)
F2y = 
figure
fplot(F2y, [-10*pi 10*pi])
grid
F2z = int(wilberforce_z * exp(1j*omega*t), t, -T, T)
F2z = 
F2z = vpa(simplify(F2z, 500), 5)
F2z = 
F2z = subs(F2z, T, 1)
F2z = 
figure
fplot(F2z, [-10*pi 10*pi])
grid
produces a more tractible result, without the δ functions.
.
  3 Comments
Luke Weston
Luke Weston on 3 Nov 2021
oops, forgot to initialise the variable... ignore :^)
Star Strider
Star Strider on 3 Nov 2021
As always, my pleasure!
(It’s necessary to add ‘omega’ and ‘T’ to the original syms call.)

Sign in to comment.

More Answers (2)

Paul
Paul on 3 Nov 2021
Replicating the code:
clear
syms t
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
k = 4.5682;
m = 0.2;
w = 0.0014;
J = 0.0327;
I = 0.00362;
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
Check the solutions, using vpa to better see their form:
vpa(qwe.y,5)
ans = 
vpa(qwe.z,5)
ans = 
So both soutions are simple sums of two cos functions. The Fourier transform of a cos() is the sum of two Dirac delta functions. So we should expect the Fourier transform of each solution to be the sum of four Dirac delta functions
vpa(fourier(qwe.y),5)
ans = 
vpa(fourier(qwe.z),5)
ans = 
which is exactly what they are.
What additional analysis might you be interested in?

Walter Roberson
Walter Roberson on 3 Nov 2021
clear
syms t
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
k = 4.5682;
m = 0.2;
w = 0.0014;
J = 0.0327;
I = 0.00362;
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
wilberforce_y = qwe.y
wilberforce_y = 
wyf = fourier(wilberforce_y)
wyf = 
wilberforce_z = qwe.z
wilberforce_z = 
wfz = fourier(wilberforce_z)
wfz = 
fplot(abs(wyf))
fplot(abs(wfz))

Community Treasure Hunt

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

Start Hunting!