Numerical integration of improper integral

[nkinar]

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.

<br /> \[<br /> 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\}<br /> \]<br />


In the equation above, {\mathop{\rm Re}\nolimits} indicates the real part, \pi is the ubiquitous pi constant, i denotes an imaginary number, and \omega is the angular frequency (1/s).

Also,

<br /> \[<br /> S\left( \omega \right) = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}\exp \left[ { - \frac{{\omega ^2 }}{{\omega _0^2 }}} \right]<br /> \]<br />

In the equation above, \omega_0 is a reference angular frequency (1/s).

Also,

<br /> \[<br /> kr = \left( {1 - \frac{i}{{2Q}}} \right)\left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \omega t<br /> \]<br />

In the equation above, Q is a constant real number, \[\gamma = (\pi Q)^{ - 1} \], and t is the time (s).

How would I numerically integrate the improper integral to obtain u(t), also using the other formulas listed above?

 

[CompuChip]

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

u(t) = Re\left\{ \int_0^\infty f(\omega) \, d\omega \right\}
where the integrand looks like
f(\omega) = C \exp( \alpha(\omega) ) \exp( i \beta(\omega))
for some real functions alpha and beta. Then you are interested in
u(t) = C \int_0^\infty e^{\alpha(\omega)} \cos(\beta(\omega)) \, d\omega
If you choose values for all your constants (Q, w₀, 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?)

 

[nkinar]

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:

<br /> \[<br /> 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\}<br /> \]<br />

Now this is close to the form that you give, but with an additional K(\omega) term:

<br /> \[<br /> f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right) - K\left( \omega \right)} \right]<br /> \]<br />

Where:

<br /> \[<br /> C = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}<br /> \]<br />


<br /> \[<br /> \alpha \left( \omega \right) = \frac{{ - \omega ^2 }}{{\omega _0^2 }}<br /> \]<br /> <br />


<br /> \[<br /> \beta \left( \omega \right) = \omega t - \omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } <br /> \]<br />


<br /> \[<br /> K\left( \omega \right) = \frac{1}{{2Q}}\omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } <br /> \]<br /> <br />


However, I would like to obtain the following:

<br /> \[<br /> f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right]<br /> \]<br />

Now using Euler's formula:

<br /> \[<br /> 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\}<br /> \]<br />


<br /> \[<br /> f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right]<br /> \]<br />

Now how do I get rid of the K(\omega) term?

 

[nkinar]

Ah - use identities!

<br /> \[<br /> 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]<br /> \]<br />

 

[nkinar]

Now the real part of the RHS becomes:

<br /> <br /> \[<br /> 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]<br /> \]

 

[CompuChip]

Yes
<br /> e^{a} e^{ib + c} = e^a e^{ib} e^{c} = e^{a + c} e^{i b}<br />

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

 

[nkinar]

That's terrific, thanks CompuChip!

Sure, normally

<br /> 0 \leq \tau \leq 2<br />

<br /> 0 \le Q \le 500<br />

<br /> \omega_0 = 2 \pi f_0<br />

<br /> 0 \leq f_0 \leq 1000<br />

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 \omega\rightarrow \infty, and become extremely large as \omega\rightarrow 0. I suspect a small error in my code. What does your Mathematica notebook say for \tau=0.2, Q=400, and \omega_0= 2\pi(50)?

 

[nkinar]

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 S(\omega) can be defined as:

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

Now by taking the FFT of zero-padded discrete s(t) and shifting the zero-frequency component to the center of the spectrum, we obtain S(\omega).

But then again, it's always conceptually nicer to do this with S(\omega) rather than go through the hassle of the FFT.

 

[nkinar]

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 K(\omega) term, which I write again below:

<br /> K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}}<br />

Now

<br /> \lim_{\omega \to \infty}K(\omega) \rightarrow \infty<br />

since \omega^{1 - 1/(\pi Q)} \rightarrow \omega as Q \rightarrow \infty, with \omega_0 \in\Re and \tau \in \Re as constants.

But then

<br /> \mbox{exp} \left( -K(\omega) \right) \rightarrow 0<br />

as \omega \rightarrow \infty.

Strangely enough, I've found that by removing the K(\omega) 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:

<br /> \[<br /> 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]<br /> \]<br />

Here is the integrand after removing the exponential term with K(\omega):

<br /> \[<br /> f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right]<br /> \]<br /> <br />

One lingering question remains: would it be logical to suppose that the contribution of \mbox{exp} \left[ -K(\omega) \right] to the numerical integral is minimal?

 

[nkinar]

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

I was wrong. After checking the code over, I traced the problem back to the following line of Matlab code used to calculate K(\omega):

<br /> K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}}<br />

This is the wrong line of Matlab code:

K = ( 1.0 / 2 * Q ) .* omega .* tau .* term;

This is the proper line:

K = ( 1.0 / ( 2 * Q )  ) .* omega .* tau .* term;

I've found that \mbox{exp}\left[ -K(\omega) \right ] is indeed required in the integrand since this seems to govern the attenuation of the wavelet.

For the integrand f(\omega), I've found that a discontinuity occurs near \omega = 0, so I had to numerically integrate using trapezoidal integration from \omega = 2\pi(0.1) to \omega = 2\pi(1000), instead of on the interval \left [ 0 , \infty \right)

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!