Evaluating integral in manual inverse fourier transform

[Tom van Woudenberg]

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(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:

Kind regards,

Tom van Woudenberg

Delft University of Technology

[Alan Bromborsky]

See attached ft.html (custom jupyter notebook image, Format and gprint not in sympy distribution yet) you should be able to calculate your integral (inverse fourier transform) analytically with the residue theorem -

https://en.wikipedia.org/wiki/Residue_theorem

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

[Alan Bromborsky]

Also note if you use the residue theorem there is a different path depending on whether t>0 or t<0.  This means when you calculate the residue you only use one pole (phi+ or phi-) depending on the path.  The path chosen makes sure that the part of it not on the real axis does not contribute to the integral.

On 2/6/23 10:36 AM, 'Tom van Woudenberg' via sympy wrote:

[Alan Bromborsky]

You can use residues but I made a mistake.  You have to do a partial fraction decomposition of the reciprocal quadratic factor.  Then both poles of your integrand lie in the lower half plane so that the inverse transform is zero for t<0.  For t>0 you have calculate the residue of each component of the partial fraction.

On 2/6/23 10:36 AM, 'Tom van Woudenberg' via sympy wrote:

[Alan Bromborsky]

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:

[Tom van Woudenberg]

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:

[Alan Bromborsky]

The best way is shown in the attached document (work in progress not done yet).  I am doing thing because I have a fascination with Fourier transforms.  You case is in Sympy.pdf.  Fourier.pdf is a crib sheet I wrote to help a friend with the Fourier transforms of modulated waveforms.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/e27ffb73-aad8-4bb1-a004-fbe0a27b9074n%40googlegroups.com.

[Alan Bromborsky]

Attached is the latest version.  This should do it.  Any questions welcomed.  Note that one of the major problems with taking inverse Fourier transforms is that in may cases (yours is an example) you wind up with piecewise functions.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/e27ffb73-aad8-4bb1-a004-fbe0a27b9074n%40googlegroups.com.

[Alan Bromborsky]

I didn't proof read well enough.  Typo in equation 4.  Correction attached

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/e27ffb73-aad8-4bb1-a004-fbe0a27b9074n%40googlegroups.com.

[Tom van Woudenberg]

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

Do you have any ideas on the difference?

Op woensdag 8 februari 2023 om 01:10:05 UTC+1 schreef brombo:

[Alan Bromborsky]

Is the plot shown here a numerical integration of your entire function or is it a plot of my H(t)?  I am going to redo the convolution using a package for convolution of piecewise functions I wrote for sympy.  I love doing this sort of stuff in my retirement.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/723c7d86-1c6d-492a-9f86-16978b4b837bn%40googlegroups.com.

[Alan Bromborsky]

I think what I have done is correct up to the point of taking the convolution product.  Both G(t) and H(t) are piecewise functions.  You have to be super careful when convolving them.  I think that is where I messed up.  I am working on a sympy class to convolve piecewise functions.  Currently it works for some but I hard coded some things in the class I shouldn't have and am currently fixing that problem.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/723c7d86-1c6d-492a-9f86-16978b4b837bn%40googlegroups.com.

[Alan Bromborsky]

Attached are latest results (I had calculated the roots of the quadratic wrong) and a plot -

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/723c7d86-1c6d-492a-9f86-16978b4b837bn%40googlegroups.com.

[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:

[Tom van Woudenberg]

This is the result in Python (same as in Maple): 

Op vrijdag 10 februari 2023 om 17:31:50 UTC+1 schreef Tom van Woudenberg:

[Alan Bromborsky]

I will make sure I translated you code to the correct fromula.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2dc42a4a-83eb-451e-b432-ef3146e076f6n%40googlegroups.com.

[Alan Bromborsky]

In your numerical solution your integration statement (I am not familiar with Maple or is this sympy where "import sympy as sp" is used) is -

solution_numeric = 1 / sp.pi * sp.integrate(sp.re(solution_in_frequency_domain_numeric*sp.exp(sp.I*2*phi*t)),(phi,0,4))

Does that statement mean you are integrating between 0 and 4 and not -infinity to infinity?  If so why did you pick the interval [0,4] to integrate?

Also why only use the real part?  A good check of the calculation is if the imaginary part is zero.

[Alan Bromborsky]

I corrected a mistake in calculating the roots of the quadratic which of course changes the time behavior of the waveform.  My numerical Fourier transform agrees with yours and I trust it.  I have to track down a multiplicative factor of pi.  The last figure shows the numerical result in blue and the analytic result times pi in red.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2dc42a4a-83eb-451e-b432-ef3146e076f6n%40googlegroups.com.

[Alan Bromborsky]

I have tracked down the errors and the analytic and numerical results now agree (see attached) -

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2dc42a4a-83eb-451e-b432-ef3146e076f6n%40googlegroups.com.

[Tom van Woudenberg]

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:

[Alan Bromborsky]

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

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/0d10619b-0d04-4101-858d-9e40d4ae1b59n%40googlegroups.com.

[Alan Bromborsky]

I am attaching a general crib sheet for Fourier Transforms much about calculating the Fourier Transform correctly from the Discreet Fourier Transform (which can be evaluated with an FFT).  Also what can be done partially analytically if you have a finite repetitive pulse train or are using a pulse or finite pulse train to modulate a sine waver.  I am also attaching the code for the plots if anyone is interested in that.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/0d10619b-0d04-4101-858d-9e40d4ae1b59n%40googlegroups.com.