# Numerical integration of improper integral

 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. $$$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\}$$$ 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, $$$S\left( \omega \right) = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}\exp \left[ { - \frac{{\omega ^2 }}{{\omega _0^2 }}} \right]$$$ In the equation above, $$\omega_0$$ is a reference angular frequency (1/s). Also, $$$kr = \left( {1 - \frac{i}{{2Q}}} \right)\left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \omega t$$$ 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?
 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 $$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, 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?)
 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: $$$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\}$$$ Now this is close to the form that you give, but with an additional $$K(\omega)$$ term: $$$f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right) - K\left( \omega \right)} \right]$$$ Where: $$$C = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}$$$ $$$\alpha \left( \omega \right) = \frac{{ - \omega ^2 }}{{\omega _0^2 }}$$$ $$$\beta \left( \omega \right) = \omega t - \omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma }$$$ $$$K\left( \omega \right) = \frac{1}{{2Q}}\omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma }$$$ However, I would like to obtain the following: $$$f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right]$$$ Now using Euler's formula: $$$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\}$$$ $$$f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right]$$$ Now how do I get rid of the $$K(\omega)$$ term?
 P: 75 Numerical integration of improper integral Ah - use identities! $$$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]$$$
 P: 75 Now the real part of the RHS becomes: $$$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]$$$
 Sci Advisor HW Helper P: 4,300 Yes $$e^{a} e^{ib + c} = e^a e^{ib} e^{c} = e^{a + c} e^{i b}$$ 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)$$
 P: 75 That's terrific, thanks CompuChip! Sure, normally $$0 \leq \tau \leq 2$$ $$0 \le Q \le 500$$ $$\omega_0 = 2 \pi f_0$$ $$0 \leq f_0 \leq 1000$$ 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)$$?
 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 $$S(\omega)$$ can be defined as: $$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)$$ 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.
 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 $$K(\omega)$$ term, which I write again below: $$K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}}$$ Now $$\lim_{\omega \to \infty}K(\omega) \rightarrow \infty$$ 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 $$\mbox{exp} \left( -K(\omega) \right) \rightarrow 0$$ 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: $$$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]$$$ Here is the integrand after removing the exponential term with $$K(\omega)$$: $$$f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right]$$$ 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?
P: 75
 Quote by 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)$$:

$$K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}}$$

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!

 Related Discussions Calculus & Beyond Homework 1 Calculus & Beyond Homework 7 Calculus & Beyond Homework 5 Engineering, Comp Sci, & Technology Homework 9 Calculus 6