Register to reply

Numerical integration of improper integral

Share this thread:
nkinar
#1
Mar16-10, 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?
Phys.Org News Partner Science news on Phys.org
Sapphire talk enlivens guesswork over iPhone 6
Geneticists offer clues to better rice, tomato crops
UConn makes 3-D copies of antique instrument parts
CompuChip
#2
Mar17-10, 05:31 AM
Sci Advisor
HW Helper
P: 4,300
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?)
nkinar
#3
Mar17-10, 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?

nkinar
#4
Mar17-10, 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]
nkinar
#5
Mar17-10, 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]
CompuChip
#6
Mar17-10, 02:31 PM
Sci Advisor
HW Helper
P: 4,300
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]
nkinar
#7
Mar17-10, 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]?
nkinar
#8
Mar17-10, 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 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.
nkinar
#9
Mar17-10, 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?
nkinar
#10
Mar17-10, 11:04 PM
P: 75
Quote Quote by nkinar View Post
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 [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:
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 [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!


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