Numerical integration of improper integralby nkinar Tags: improper, integral, integration, numerical 

#1
Mar1610, 11:44 PM

P: 75

Hello
I need to generate synthetic data to test an algorithm used to process data from an experiment. A synthetic wavelet is constructed using the following equations, but I am uncertain how to numerically evaluate the improper integral shown below. [tex] \[ u(t) = {\mathop{\rm Re}\nolimits} \left\{ {{\textstyle{1 \over \pi }}\int\limits_0^\infty {S\left( \omega \right)\exp \left[ {i\left( {\omega t  kr} \right)} \right]d\omega } } \right\} \] [/tex] In the equation above, [tex]{\mathop{\rm Re}\nolimits} [/tex] indicates the real part, [tex]\pi[/tex] is the ubiquitous pi constant, [tex]i[/tex] denotes an imaginary number, and [tex]\omega[/tex] is the angular frequency (1/s). Also, [tex] \[ S\left( \omega \right) = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}\exp \left[ {  \frac{{\omega ^2 }}{{\omega _0^2 }}} \right] \] [/tex] In the equation above, [tex]\omega_0[/tex] is a reference angular frequency (1/s). Also, [tex] \[ kr = \left( {1  \frac{i}{{2Q}}} \right)\left( {\frac{\omega }{{\omega _0 }}} \right)^{  \gamma } \omega t \] [/tex] In the equation above, [tex]Q[/tex] is a constant real number, [tex]\[\gamma = (\pi Q)^{  1} \] [/tex], and [tex]t[/tex] is the time (s). How would I numerically integrate the improper integral to obtain [tex]u(t)[/tex], also using the other formulas listed above? 



#2
Mar1710, 05:31 AM

Sci Advisor
HW Helper
P: 4,301

I took a quick look but I don't have time to really go into it.
If I were you, I would start by plugging in all the definitions, to obtain something of the form [tex]u(t) = Re\left\{ \int_0^\infty f(\omega) \, d\omega \right\} [/tex] where the integrand looks like [tex]f(\omega) = C \exp( \alpha(\omega) ) \exp( i \beta(\omega)) [/tex] for some real functions alpha and beta. Then you are interested in [tex]u(t) = C \int_0^\infty e^{\alpha(\omega)} \cos(\beta(\omega)) \, d\omega[/tex] If you choose values for all your constants (Q, w_{0}, t) you should be able to plot this as a function of omega and see how quickly it oscillates around omega = 0 and how quickly it goes to zero for omega large. Then you should choose a numerical integration strategy based on that (e.g. do you need to make approximations? Do you write a C++ program or let Mathematica's NIntegrate have a go?) 



#3
Mar1710, 12:18 PM

P: 75

Hello
Thank you very much for your response, CompuChip! Okay, here is my progress so far. By making all of the substitutions, I obtain the following for the integrand: [tex] \[ f\left( \omega \right) = \left\{ {4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}\exp \left[ {\frac{{  \omega ^2 }}{{\omega _0^2 }}} \right]\exp \left[ {i\left[ {\omega t  \omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{  \gamma } } \right]  \frac{1}{{2Q}}\omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{  \gamma } } \right]} \right\} \] [/tex] Now this is close to the form that you give, but with an additional [tex]K(\omega)[/tex] term: [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)  K\left( \omega \right)} \right] \] [/tex] Where: [tex] \[ C = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }} \] [/tex] [tex] \[ \alpha \left( \omega \right) = \frac{{  \omega ^2 }}{{\omega _0^2 }} \] [/tex] [tex] \[ \beta \left( \omega \right) = \omega t  \omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{  \gamma } \] [/tex] [tex] \[ K\left( \omega \right) = \frac{1}{{2Q}}\omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{  \gamma } \] [/tex] However, I would like to obtain the following: [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right] \] [/tex] Now using Euler's formula: [tex] \[ f\left( \omega \right) = {\mathop{\rm Re}\nolimits} \left\{ {C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right]} \right\} \] [/tex] [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \] [/tex] Now how do I get rid of the [tex]K(\omega)[/tex] term? 



#4
Mar1710, 12:24 PM

P: 75

Numerical integration of improper integral
Ah  use identities!
[tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right]\exp \left[ {  K\left( \omega \right)} \right] \] [/tex] 



#5
Mar1710, 12:26 PM

P: 75

Now the real part of the RHS becomes:
[tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {  K\left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \] [/tex] 



#6
Mar1710, 02:31 PM

Sci Advisor
HW Helper
P: 4,301

Yes
[tex] e^{a} e^{ib + c} = e^a e^{ib} e^{c} = e^{a + c} e^{i b} [/tex] Looks good. Now can you give some info about the values of the constants? If I choose them somewhat conveniently, for example [tex]\tau = 0, Q = 1/\pi, \omega_0 = 1[/tex] I can even let Mathematica do the exact integration and get [tex]u(t) = \int_0^\infty f(\omega) \, d\omega =  \frac{\sqrt{\pi}}{8} e^{t / 4} (t^2  2)[/tex] 



#7
Mar1710, 02:58 PM

P: 75

That's terrific, thanks CompuChip!
Sure, normally [tex] 0 \leq \tau \leq 2 [/tex] [tex] 0 \le Q \le 500 [/tex] [tex] \omega_0 = 2 \pi f_0 [/tex] [tex] 0 \leq f_0 \leq 1000 [/tex] That's great to know that you can do the exact integration using Mathematica! On my end here (using Matlab), I'm trying to do everything numerically, and I am finding that the values of the integrand are very small as [tex]\omega\rightarrow \infty [/tex], and become extremely large as [tex]\omega\rightarrow 0[/tex]. I suspect a small error in my code. What does your Mathematica notebook say for [tex]\tau=0.2[/tex], [tex]Q=400[/tex], and [tex]\omega_0= 2\pi(50)[/tex]? 



#8
Mar1710, 03:27 PM

P: 75

I've also been thinking that another way to go about it would be to do the following.
In the time domain, an equivalent version of [tex]S(\omega)[/tex] can be defined as: [tex] s(t) = \left(1  \frac{1}{2}\omega_0^2 t^2\right) \mbox{exp} \left(\frac{1}{4}\omega_0^2 t^2 \right) [/tex] Now by taking the FFT of zeropadded discrete [tex]s(t)[/tex] and shifting the zerofrequency component to the center of the spectrum, we obtain [tex]S(\omega)[/tex]. But then again, it's always conceptually nicer to do this with [tex]S(\omega)[/tex] rather than go through the hassle of the FFT. 



#9
Mar1710, 10:05 PM

P: 75

After checking things over a number of times, I do not think that there is any error in my Matlab code. (At least I hope not.) The difficulty seems to arise with the [tex]K(\omega)[/tex] term, which I write again below:
[tex] K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{\frac{1}{\pi Q}} [/tex] Now [tex] \lim_{\omega \to \infty}K(\omega) \rightarrow \infty [/tex] since [tex] \omega^{1  1/(\pi Q)} \rightarrow \omega [/tex] as [tex]Q \rightarrow \infty[/tex], with [tex]\omega_0 \in\Re[/tex] and [tex]\tau \in \Re [/tex] as constants. But then [tex] \mbox{exp} \left( K(\omega) \right) \rightarrow 0 [/tex] as [tex]\omega \rightarrow \infty[/tex]. Strangely enough, I've found that by removing the [tex]K(\omega)[/tex] term from the expression, I can obtain what appears to be exactly the output that I require! So here is the integrand before the removal: [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {  K\left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \] [/tex] Here is the integrand after removing the exponential term with [tex]K(\omega)[/tex]: [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \] [/tex] One lingering question remains: would it be logical to suppose that the contribution of [tex] \mbox{exp} \left[ K(\omega) \right] [/tex] to the numerical integral is minimal? 



#10
Mar1710, 11:04 PM

P: 75

[tex] K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{\frac{1}{\pi Q}} [/tex] This is the wrong line of Matlab code:
For the integrand [tex]f(\omega)[/tex], I've found that a discontinuity occurs near [tex]\omega = 0[/tex], so I had to numerically integrate using trapezoidal integration from [tex]\omega = 2\pi(0.1)[/tex] to [tex] \omega = 2\pi(1000) [/tex], instead of on the interval [tex] \left [ 0 , \infty \right) [/tex] Now finally I am able to calculate the wavelet which is required to test my algorithm! Many thanks for your help, CompuChip. I couldn't have done it without your suggestions! 


Register to reply 
Related Discussions  
Improper Integral Integration  Calculus & Beyond Homework  1  
Improper Integration  Calculus & Beyond Homework  7  
Improper integration  Calculus & Beyond Homework  5  
numerical analysis (composite numerical integration)  Engineering, Comp Sci, & Technology Homework  9  
Limits  Improper Integration  Calculus  6 