Fourier transform of symbolic function

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?

Answer 1

[Star Strider]

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.

.

[Luke Weston]

Hi, thanks for the answer, although it does seem that matlab throws an error (for me) at the line:

F2y = int(wilberforce_y * exp(1j*omega*t), t, -T, T)

and:

F2z = int(wilberforce_z * exp(1j*omega*t), t, -T, T)

the error is as follows:

"Unrecognized function or variable 'omega'."

It doesn't seem like omega is defined anywhere, what is the purpose of this variable?

[Luke Weston]

oops, forgot to initialise the variable... ignore :^)

[Star Strider]

As always, my pleasure!

(It’s necessary to add ‘omega’ and ‘T’ to the original syms call.)

Answer 2

[Paul]

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?

Answer 3

[Walter Roberson]

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