Thank you for your comment. In my write up I am proselytizing two
things. Firstly, the Asymptote software package makes really nice plots
(I will attach the code I wrote to this email). Secondly, the difference
between the FFT (inverse FFT) and the Fourier transform (inverse Fourier
transform) for numerical calculations should be made very clear. The
FFT is actually the Fourier transform of z(t) = sum(z_j*delta(t-j*dt), 0
<= j < N) where delta(t) is the delta function. The numerical Fourier
transform I calculate is the Fourier transform of z(t) =
sum(z_j*triangle_dt(t-j*dt), 0 <= j < N).
On 2/14/23 6:51 AM, 'Tom van Woudenberg' via sympy wrote:
Hi Brombo,
That's a very nice result. I'll use it in my education!
Regarding your previous questions: The correct statement in
Python/SymPy should be:
solution = 2 *
sp.integrate(sp.re(solution_in_frequency_domain*sp.exp(sp.I*2*sp.pi*phi*t)),(phi,0,4))
I arbitrarily chose 4 as an cutoff value for the frequenciesbecause
higher values seem negligible. I only used the real part because the
original problem is the real-valued vibrations of a SDOF system fored
by a block function.
Op maandag 13 februari 2023 om 16:31:29 UTC+1 schreef brombo:
I have tracked down the errors and the analytic and numerical
results now agree (see attached) -
On 2/10/23 11:32 AM, 'Tom van Woudenberg' via sympy wrote:
This is the result in Python (same as in Maple): downloaden (5).png
Op vrijdag 10 februari 2023 om 17:31:50 UTC+1 schreef Tom van
Woudenberg:
Hi Brombo,
Thank you for the update. It seems my previous posts didn't
show up. Anyway, you result doesn't match the result in Maple
and the numerical evalution of the integral in Python:
Would be wonderful if we'd find an analytical solution.
Op vrijdag 10 februari 2023 om 01:08:56 UTC+1 schreef brombo:
Attached are latest results (I had calculated the roots
of the quadratic wrong) and a plot -
On 2/8/23 4:24 AM, 'Tom van Woudenberg' via sympy wrote:
Hi Brombo,
Thank you for the extensive working-out. I really
appreciate that!
However, the result doesn't seem to match the result in
got in Maple (below result in Python for N(t):
Schermafbeelding 2023-02-08 094041.jpg
Do you have any ideas on the difference?
Op woensdag 8 februari 2023 om 01:10:05 UTC+1 schreef
brombo:
I didn't proof read well enough. Typo in equation
4. Correction attached
On 2/7/23 3:02 AM, 'Tom van Woudenberg' via sympy wrote:
Thank you brombo, I'll take a closer look at the
file you send me!
Op maandag 6 februari 2023 om 22:29:25 UTC+1
schreef brombo:
I cleaned things up here is what the notebook
looks like (see attached html) -
On 2/6/23 10:36 AM, 'Tom van Woudenberg' via
sympy wrote:
Hi there,
When trying to solve a integral as part of a
manual inverse fourier transform, SymPy return
the unevaluated integral. Does anybody know if
SymPy is able to solve this integral with some
help? It would be good enough if I'd be able
to obtain the result for specific values of t.
import sympy as sp
phi,t = sp.symbols('phi,t',real=True)
sp.I*(1 -
sp.exp(4*sp.I*sp.pi*phi))*sp.exp(-8*sp.I*sp.pi*phi)/(2*sp.pi*phi*(-4*sp.pi**2*phi**2
+ 1.5*sp.I*sp.pi*phi + 4))
solution_numeric = 1 / sp.pi *
sp.integrate(sp.re
<http://sp.re>(solution_in_frequency_domain_numeric*sp.exp(sp.I*2*phi*t)),(phi,0,4))
print(solution_numeric)
returns:
(Integral(sin(4*pi*phi)*re(exp(2*I*phi*t)/(-4*pi**2*phi**2*exp(8*I*pi*phi)
+ 1.5*I*pi*phi*exp(8*I*pi*phi) +
4*exp(8*I*pi*phi))), (phi, 0, 4)) +
Integral(cos(4*pi*phi)*im(exp(2*I*phi*t)/(-4*pi**2*phi**2*exp(8*I*pi*phi)
+ 1.5*I*pi*phi*exp(8*I*pi*phi) +
4*exp(8*I*pi*phi))), (phi, 0, 4)) +
Integral(-im(exp(2*I*phi*t)/(-4*pi**2*phi**2*exp(8*I*pi*phi)
+ 1.5*I*pi*phi*exp(8*I*pi*phi) +
4*exp(8*I*pi*phi))), (phi, 0, 4)))/(2*pi**2*phi)
Plotting the result for t,0,15 should give
this result according to Maple:
Schermafbeelding 2023-02-06 163521.jpg
Kind regards,
Tom van Woudenberg
Delft University of Technology
--
You received this message because you are
subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop
receiving emails from it, send an email to
sympy+un...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/eea7eaef-8752-41f8-bf9d-ba78a1782c37n%40googlegroups.com
<https://groups.google.com/d/msgid/sympy/eea7eaef-8752-41f8-bf9d-ba78a1782c37n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
You received this message because you are
subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving
emails from it, send an email to
sympy+un...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/e27ffb73-aad8-4bb1-a004-fbe0a27b9074n%40googlegroups.com
<https://groups.google.com/d/msgid/sympy/e27ffb73-aad8-4bb1-a004-fbe0a27b9074n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
You received this message because you are subscribed to
the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails
from it, send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/723c7d86-1c6d-492a-9f86-16978b4b837bn%40googlegroups.com
<https://groups.google.com/d/msgid/sympy/723c7d86-1c6d-492a-9f86-16978b4b837bn%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
You received this message because you are subscribed to the
Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it,
send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/2dc42a4a-83eb-451e-b432-ef3146e076f6n%40googlegroups.com
<https://groups.google.com/d/msgid/sympy/2dc42a4a-83eb-451e-b432-ef3146e076f6n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
You received this message because you are subscribed to the Google
Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send
an email to sympy+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/0d10619b-0d04-4101-858d-9e40d4ae1b59n%40googlegroups.com
<https://groups.google.com/d/msgid/sympy/0d10619b-0d04-4101-858d-9e40d4ae1b59n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sympy+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/324f9338-193e-07e3-d381-1fb15aeffd02%40gmail.com.
// To use this code you need to install the fftw library
import graph;
import math;
struct real_f_data
{
real [] t;
real [] f;
}
pair a,b,c;
a = (1,0);
b = (0,-3/4);
c = (-4,0);
pair wp = (-b+sqrt(b**2-4*a*c))/(2*a);
pair wm = (-b-sqrt(b**2-4*a*c))/(2*a);
//real wi = wp.y;
//real wr = wp.x;
real wi = 3/8;
real wr = sqrt(247)/8;
real wabs = sqrt(wi**2+wr**2);
real wth = atan2(wi,wr);
real dw = 2*wr;
real dw_inv = 1/(2*wr);
real sinc(real x)
{
if ( x != 0.0 )
{
return sin(x)/x;
}
return 1.0;
}
real P(real t)
{
return -2*exp(-wi*t)*cos(wr*t-wth)/(wabs*dw);
}
real Nanalytic(real t)
{
if ((t >= 2.0) && (t <= 4.0))
{
return (P(t-2)-P(0));
}
if (t > 4.0)
{
return (P(t-2)-P(t-4));
}
return 0.0;
}
pair gateFT(real w)
{
real x,y;
if (w != 0.0)
{
x = 2*sin(w)/w;
y = 0.0;
}
else
{
x = 2.0;
y = 0.0;
}
return (x,y);
}
pair [] F_pair_array(pair F(real w), real dw, int N)
{
pair [] Z;
for (int i = 0; i < N; ++i)
{
Z.push(F(i*dw));
}
return Z;
}
real Re_gateFT(real w)
{
return gateFT(w).x;
}
real_f_data Inverse_Fourier_Transform(pair [] F,real dw)
{
int N = F.length-1;
real dt = 2*pi/(N*dw);
real pi_div_N = pi/N;
pair [] f_INVFFT = fft(F,sign=1);
real [] f,t;
real tmp;
for (int i = 0; i<= N; ++i)
{
tmp = pi_div_N*i;
t.push(i*dt);
f.push(dw*sinc(tmp)**2*f_INVFFT[i].x/pi);
}
real_f_data f_t;
f_t.t = t;
f_t.f = f;
return f_t;
}
pair N_of_w(real w)
{
return -2*sinc(w)*exp(-3*I*w)/((w**2,0)-3*I*w/4-(4,0));
}
picture pic1;
size(pic1,250,200,IgnoreAspect);
draw(pic1,graph(Nanalytic,0,25,5000),blue,"$N(t)$");
yequals(pic1,0.0,black);
xlimits(pic1,0,25);
ylimits(pic1,-0.2,0.5);
xaxis(pic1,"$t$",Bottom,LeftTicks(step=1.0));
yaxis(pic1,"$\hat{N}(t)$",LeftRight,RightTicks(trailingzero));
shipout("pic1",pic1);
picture pic2;
size(pic2,250,200,IgnoreAspect);
draw(pic2,graph(Re_gateFT,-30*pi,30*pi,5000),blue,"$F(\omega)$");
yequals(pic2,0.0,black);
xaxis(pic2,"$\omega$",Bottom,LeftTicks);
yaxis(pic2,"$F(\omega)$",LeftRight,RightTicks(trailingzero));
shipout("pic2",pic2);
int N = 2**14;
real dw = 60*pi/N;
pair [] Z = F_pair_array(gateFT,dw,N);
real_f_data f_of_t = Inverse_Fourier_Transform(Z,dw);
picture pic3;
size(pic3,250,200,IgnoreAspect);
draw(pic3,graph(f_of_t.t[0:150],f_of_t.f[0:150]),blue,"$f(t)$");
xaxis(pic3,"$t$",Bottom,LeftTicks);
yaxis(pic3,"$f(t)$",LeftRight,RightTicks(trailingzero));
shipout("pic3",pic3);
int N = 2**14;
real dw = 60*pi/N;
pair [] M = F_pair_array(N_of_w,dw,N);
real_f_data N_of_t = Inverse_Fourier_Transform(M,dw);
picture pic4;
size(pic4,250,200,IgnoreAspect);
draw(pic4,graph(N_of_t.t[0:750],N_of_t.f[0:750]),blue+linewidth(2),"Numeric");
draw(pic4,graph(Nanalytic,0,25,5000),red,"Analytic");
yequals(pic4,0.0,black);
xaxis(pic4,"$t$",Bottom,LeftTicks(step=1));
yaxis(pic4,"$\hat{N}(t)$",LeftRight,RightTicks(trailingzero),ymin=-0.2,ymax=0.4);
add(pic4,legend(pic4),point(pic4,NE),25SW,UnFill);
shipout("pic4",pic4);
//More tests
pair G_of_w(real w)
{
return -2*sinc(w)*exp(-3*I*w);
}
int N = 2**14;
real dw = 60*pi/N;
pair [] Gw = F_pair_array(G_of_w,dw,N);
real_f_data G_of_t = Inverse_Fourier_Transform(Gw,dw);
picture pic5;
size(pic5,250,200,IgnoreAspect);
draw(pic5,graph(G_of_t.t[0:500],G_of_t.f[0:500]),blue,"$G(t)$");
xaxis(pic5,"$t$",Bottom,LeftTicks(step=1));
yaxis(pic5,"$\hat{G}(t)$",LeftRight,RightTicks(trailingzero));
shipout("pic5",pic5);
pair H_of_w(real w)
{
return 1/((w**2,0)-3*I*w/4-(4,0));
}
int N = 2**14;
real dw = 60*pi/N;
pair [] Hw = F_pair_array(H_of_w,dw,N);
real_f_data H_of_t = Inverse_Fourier_Transform(Hw,dw);
picture pic6;
size(pic6,250,200,IgnoreAspect);
draw(pic6,graph(H_of_t.t[0:500],H_of_t.f[0:500]),blue,"$H(t)$");
xaxis(pic6,"$t$",Bottom,LeftTicks(step=1));
yaxis(pic6,"$\hat{H}(t)$",LeftRight,RightTicks(trailingzero));
shipout("pic6",pic6);
pair H1_of_w(real w)
{
return dw_inv*(1/(w-wp)-1/(w-wm));
}
int N = 2**14;
real dw = 60*pi/N;
pair [] Hw1 = F_pair_array(H1_of_w,dw,N);
real_f_data H1_of_t = Inverse_Fourier_Transform(Hw1,dw);
picture pic7;
size(pic7,250,200,IgnoreAspect);
draw(pic7,graph(H1_of_t.t[0:500],H1_of_t.f[0:500]),blue,"$H(t)$");
xaxis(pic7,"$t$",Bottom,LeftTicks(step=1));
yaxis(pic7,"$\hat{H}(t)$",LeftRight,RightTicks(trailingzero));
shipout("pic7",pic7);
