Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerical integration of improper integral

  1. Mar 16, 2010 #1
    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. jcsd
  3. Mar 17, 2010 #2

    CompuChip

    User Avatar
    Science Advisor
    Homework Helper

    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, w0, 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?)
     
  4. Mar 17, 2010 #3
    Hello--

    Thank you very much for your response, CompuChip! :smile:

    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?
     
  5. Mar 17, 2010 #4
    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]
     
  6. Mar 17, 2010 #5
    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]
     
  7. Mar 17, 2010 #6

    CompuChip

    User Avatar
    Science Advisor
    Homework Helper

    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]
     
  8. Mar 17, 2010 #7
    That's terrific, thanks CompuChip! :smile:

    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. :eek: What does your Mathematica notebook say for [tex]\tau=0.2[/tex], [tex]Q=400[/tex], and [tex]\omega_0= 2\pi(50)[/tex]?
     
  9. Mar 17, 2010 #8
    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 zero-padded discrete [tex]s(t)[/tex] and shifting the zero-frequency 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.
     
  10. Mar 17, 2010 #9
    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? :confused:
     
  11. Mar 17, 2010 #10
    I was wrong. :redface: After checking the code over, I traced the problem back to the following line of Matlab code used to calculate [tex]K(\omega)[/tex]:

    [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:
    Code (Text):

    K = ( 1.0 / 2 * Q ) .* omega .* tau .* term;
     
    This is the proper line:

    Code (Text):

    K = ( 1.0 / ( 2 * Q )  ) .* omega .* tau .* term;
     
    I've found that [tex]\mbox{exp}\left[ -K(\omega) \right ] [/tex] is indeed required in the integrand since this seems to govern the attenuation of the wavelet.

    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! :smile: :smile: :smile:
     
    Last edited: Mar 17, 2010
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Similar Discussions: Numerical integration of improper integral
  1. Improper Integration (Replies: 5)

  2. Improper integrals (Replies: 4)

  3. Improper Integrals (Replies: 1)

  4. "Improper" Integral (Replies: 3)

Loading...