Integrating a product of exponential and trigonometric functions

[Talon44]

TL;DR Summary

Trying to solve an integral involving an exponential function and trigonometric function. Looking for some tips.

I am looking for a closed form solution to an integral of the form:

$$ \int_0^\infty \frac{e^{-Du^2t}u \sin{ux}}{u^2+h^2} du $$

D, t, and h are positive and x is unrestricted.

I have tried everything, integration by parts, substitution, even complex integration with residue analysis. I've gotten glimmers of something intelligible but always end up at a dead end. I'm looking for some help seeing some strategy I might be missing.

Happy to provide some more detail about what I've tried but honestly I've got so many pages of scribbles here I'm not even sure what I would include at this point.

 

[anuttarasammyak]

Residue theorem in complex integral seems to provide you the result.

 

[Talon44]

Thank you for your response.

Actually, what I'm really trying to do is show that

$$ \frac{2hD}{\pi} \int_0^\infty \frac{e^{-Du^2t}u \sin{ux}}{u^2+h^2} du + \frac{2D}{\pi} \int_0^\infty \frac{e^{-Du^2t}u^2 \cos{ux}}{u^2+h^2} du $$

Has the analytical solution

$$ \sqrt{\frac{D}{\pi t}} e^{-x^2/4Dt}- hD e^{hx+Dth^2} erfc(\frac{x}{2 \sqrt{Dt}}+h\sqrt{Dt}) $$

The latter expression comes from a table of Laplace transforms. As a matter of practice for more difficult examples which do not appear in the table, I am trying to derive it myself. I integrated numerically for a few cases using Mathematica so the expressions do appear to be equivalent; I am just struggling to show it.

I did try using the residue theorem using a semicircular contour that extends from negative infinity to infinity along the real axies, then circles around the positive imaginary axis, enclosing one of the single poles at z = ih. As the radius of this semicircle extends to infinity, the integral along the real axis from negative infinity to positive infinity reduces to twice the integral that I'm looking for (which I believe is an even function). I assumed (possibly wrongly, but I didn't know how to evaluate it in any case) that the integral along the circular contour goes to zero as R goes to infinity, leaving only the residue at the single pole. But that didn't give me anything remotely like the solution. So clearly I'm doing something wrong.

 

[anuttarasammyak]

I assumed (possibly wrongly, but I didn't know how to evaluate it in any case) that the integral along the circular contour goes to zero as R goes to infinity, leaving only the residue at the single pole. But that didn't give me anything remotely like the solution. So clearly I'm doing something wrong.

Jordan’s lemma assures it.

 

[Talon44]

Ok, I'll play along then. Here is what I did.

For

$$ \frac{2hD}{\pi} \int_0^\infty \frac{e^{-Du^2t}u \sin{ux}}{u^2+h^2} du $$

I substituted in $e^{iux}$ for the sine function, expecting only to use the imaginary part of the solution, so

$$ f(z)= \frac{e^{-Dz^2t} e^{izx}z}{(z+ih)(z-ih)} $$

So then with the circular contour going to zero at R approaching infinity, I am left with

$$ \int_{-\infty}^{\infty} f(z) dz = 2 \int_0^{\infty} f(z) dz = 2 \pi i Res_{z=ih}(f(z)) $$

$$ Res_{z=ih}(f(z)) = \frac{e^{-D(ih)^2t+i(ih)x} (ih)}{2ih} = \frac{e^{Dh^2t-hx}}{2} $$

Therefore

$$ \int_0^{\infty} f(z) dz = \frac {\pi i}{2} e^{Dh^2t-hx} $$

My integral of interest is only the imaginary component, so I arrive at

$$ \frac{2hD}{\pi} \int_0^\infty \frac{e^{-Du^2t}u \sin{ux}}{u^2+h^2} du = hD e^{Dh^2t-hx} $$

I get pretty much the same result with the cosine integral, which itself feels wrong to me. So, I feel I must be doing something fundamentally incorrect. Never mind the fact that it isn't even close to what the table of transforms tells me is the solution. I am no expert at complex analysis.

 

[anuttarasammyak]

I substituted in eiux for the sine function, expecting only to use the imaginary part of the solution, so

Could you let me know whether the substitution
$$\sin ux=\frac{e^{iux}-e^{-iux}}{2i}$$
gives the same result ? For cos
$$\cos ux=\frac{e^{iux}+e^{-iux}}{2}$$

 

[Talon44]

Just to be clear, the substitutions I originally made were:

$$ \sin(ux) = \Im(e^{iux}) $$
$$ \cos(ux) = \Re(e^{iux}) $$

However, I tried it with the alternative substitutions you recommended and while the intermediate solutions to the residues for the sin and cos integrals were different, the overall result for the sum of the two integrals was not affected by the type of substitution. (Well, I've lost a factor of 2 somewhere, but the functional forms are identical.) I certainly am not getting the ERFC factor or the leading term with $e^{-x^2/4Dt}$.

 

[pasmith]

It looks like your first integral arises from a solution of the heat equation. But what initial and boundary conditions are you using?

 

[renormalize]

The latter expression comes from a table of Laplace transforms.

Can you cite the source for this table of Laplace transforms?

 

[Talon44]

J Crank, The Mathematics of Diffusion 2nd Ed. Table 2.2. A similar transform (although with heat transfer variables) appears in Carslaw and Jaeger, Conduction of Heat in Solids.

The transform I am practicing with is this:
If
$$ \bar{v}(p) = \frac{e^{-qx}}{q+h} $$
Then
$$ v(t) = \sqrt{\frac{D}{\pi t}} e^{-x^2/4Dt}- hD e^{hx+Dth^2} erfc\left(\frac{x}{2 \sqrt{Dt}}-h\sqrt{Dt}\right) $$

$p$ is the transform variable, $D$ is the diffusion coefficient, $t$ is time, $x$ is position, and $h$ is a (positive) constant. $q = \sqrt{p/D}$.

To perform the inverse transform, I used a Bromwich keyhole contour, a method outlined in Carslaw and Jaeger to inverse transform solutions to the subsidiary equation that have a branch cut. This is how I got to the integral expression that has me stuck:

$$ v(t)=\frac{2hD}{\pi} \int_0^\infty \frac{e^{-Du^2t}u \sin{ux}}{u^2+h^2} du + \frac{2D}{\pi} \int_0^\infty \frac{e^{-Du^2t}u^2 \cos{ux}}{u^2+h^2} $$

I’ll note that I’ve used this method successfully for some of the other transforms in the Crank’s table, so I’m pretty confident of my method. This one was just given me an integral I cannot seem to solve, even though it does have an analytical solution. Also as noted earlier numeric evaluation of the integral expression with a few test values of h, D, x and t does appear to return the transformed expression, so I believe my integral expression is correct. I just am struggling to show it for general values of h, D, x and t.

@pasmith

I'm not trying to solve a particular diffusion problem so I don't have any specific boundary conditions. I just trying to reproduce the inverse transform listed in Crank's table. However, my recollection is that cases like can arise when diffusion occurs through a thin boundary layer, which provides a resistance to the diffusion flux at the boundary, in which case h would be a proportionality constant.

 

[renormalize]

The transform I am practicing with is this:
If
$$ \bar{v}(p) = \frac{e^{-qx}}{q+h} $$
Then
$$ v(t) = \sqrt{\frac{D}{\pi t}} e^{-x^2/4Dt}- hD e^{hx+Dth^2} erfc\left(\frac{x}{2 \sqrt{Dt}}-h\sqrt{Dt}\right) $$

$p$ is the transform variable, $D$ is the diffusion coefficient, $t$ is time, $x$ is position, and $h$ is an unrestricted constant. $q = \sqrt{p/D}$.

I don't know that I can offer much help since I'm not particularly practiced with Laplace transforms. That said, it would be helpful to better understand your question if you could write-out the transform-pairs that you're trying to verify, explicitly as equations involving integrals that clearly display the integrand, integration variable and integral limits and/or integration contour.

 

[Talon44]

I appreciate your reply. At this point the transform is really just context. I am operating under the assumption that I have (inverse) transformed the solution to the subsidiary equation correctly. The integral form I am stuck at is already a transformed expression, one that is perfectly usable provided I am willing to do a numerical integration. However as a practical matter the analytical form in the transform table is much easier to use with data fitting programs like Excel and Origin because generally we are trying to determine a diffusion coefficient from concentration v time data.

Anyway, that's all to say that my problem isn't with the transform - it's with converting the integral representation of it that I derived into the analytical one that's listed in a table. I have an integral (well, the sum of two integrals) and I know what the analytical solution is (from the transform table) but I do not know how to get from the one to the other using standard tools of integration. This is most explicitly laid out in my second post - I'm not sure how to put my question more directly than that.

(I realize the exercise is academic. But I have some other problems I need analytical solutions to that do not appear in transform tables, and their integral representations are similar, though more complex, than this one.)

 

[anuttarasammyak]

My integral of interest is only the imaginary component, so I arrive at

2hDπ∫0∞e−Du2tusin⁡uxu2+h2du=hDeDh2t−hx

It is obvious that for x=0, the integral is zero. Your result does not satisfy it.

 

[Talon44]

Well, yes, it was already obvious that my result was not correct, which is why I am here. I know what the correct answer is - I'm trying to figure out how to derive it. It would be more useful if you could tell me what I'm doing wrong.

 

[pasmith]

Start with a simpler problem: can you show that the inverse laplace transform of $$\frac{1}{\sqrt{p} + h}$$ is that given by Wolfram Alpha, $$\frac{1}{\sqrt{\pi t}} - h e^{h^2 t} \operatorname{erfc}(h \sqrt{t})$$ Then you have $$\frac{1}{\sqrt{p/D} + h} = \frac{\sqrt{D}}{\sqrt{p} + \sqrt{D}h}.$$

 

[Talon44]

Yes, it is straightforward to use the residue theorem with the Bromwich line integral to show that if

$$ \bar{v}(p) = \frac{1}{\sqrt{p}+h} $$

Then

$$ v(t) = \frac{2}{\pi} \int_0^{\infty} {\frac{e^{-tu^2} u^2 du}{h^2+u^2}} $$

Which has the analytical solution (according to Mathematica)

$$ v(t) = \sqrt{\frac{1}{\pi t}} - h e^{h^2t} erfc(h\sqrt{t}) $$

(Subject to certain conditions, like h, t being positive.)

But as I noted, performing the inverse transform is not my problem. It’s getting from the integral representation to the closed form representation that is my challenge. This integral is simpler; Mathematica can do it. But Mathematica is unable to give me an analytical solution for the integral representation in my more difficult problem… which means I have to do it manually. That’s what I can't figure out how to show.

 

[pasmith]

If you can't figure out how to do it yourself - without relying on mathematica - in this simpler case, how do you expect to be able to adapt those steps to a case that mathematica can't handle?

The result we are looking for is $$\frac{2}{\pi} \int_0^\infty e^{-u^2 t}\frac{u^2}{u^2 + h^2}\,du = \frac1{\sqrt{\pi t}} - \frac{2he^{h^2t}}{\sqrt{\pi}} \int_{h\sqrt{t}}^\infty e^{-v^2} \,dv.$$

Substituting $v^2 = t(h^2 + u^2)$ gives $$\begin{split} \frac{2}{\pi} \int_0^\infty e^{-u^2 t}\frac{u^2}{u^2 + h^2}\,du &= \frac{2e^{h^2t}}{\pi \sqrt{t}}\int_{h\sqrt{t}}^\infty e^{-v^2} \sqrt{1 - \frac{h^2 t}{v^2}}\,dv\\ &= \frac{2he^{h^2t}}{\pi} \int_{h\sqrt{t}}^\infty e^{-v^2}\frac 1v \sqrt{ \frac{v^2}{h^2t} - 1}\,dv\end{split}$$ and now integration by parts can come into play.

 

[renormalize]

@Talon44, it turns out that we can use Mathematica to perform your integrations by borrowing a page from Feynman's book of clever integration tricks and introducing an additional parameter under the integral, followed by differentiating with respect to it. Begin by rewriting your sum of two integrals as:$$I_{\text{Talon44}}\equiv\frac{2hD}{\pi}\int_{0}^{\infty}\frac{e^{-Dtu^{2}}u\sin\left(ux\right)}{u^{2}+h^{2}}du+\frac{2D}{\pi}\int_{0}^{\infty}\frac{e^{-Dtu^{2}}u^{2}\cos\left(ux\right)}{u^{2}+h^{2}}du$$$$=-\frac{2D}{\pi}\left(h\frac{\partial}{\partial x}+\frac{\partial^{2}}{\partial x^{2}}\right)I_{\mathrm{M}}\tag{1}$$in terms of a single "master integral" $I_{\mathrm{M}}$:$$I_{\mathrm{M}}\equiv\int_{0}^{\infty}\frac{e^{-\alpha u^{2}}\cos\left(ux\right)}{u^{2}+h^{2}}du=e^{\alpha h^{2}}\int_{0}^{\infty}\frac{e^{-\alpha\left(u^{2}+h^{2}\right)}\cos\left(ux\right)}{u^{2}+h^{2}}du,\;\alpha\equiv Dt\tag{2a,b}$$We generalize this expression by introducing a Feynman parameter $p$:$$I_{\mathrm{M}}\left(p\right)\equiv e^{\alpha h^{2}}\int_{0}^{\infty}\frac{e^{-p\alpha\left(u^{2}+h^{2}\right)}\cos\left(ux\right)}{u^{2}+h^{2}}du\tag{3}$$We're of course interested in the particular value $I_{\mathrm{M}}(1)$ at $p=1$, but we'll also need the fact that $I_{\mathrm{M}}(\infty)=0$ (this is because the integrand in (3) vanishes as $p\rightarrow\infty$, so long as $\alpha\equiv Dt>0$). Differentiating eq.(3) with respect to $p$ yields an integral that can be carried out by Mathematica:$$\frac{d I_{\mathrm{M}}\left(p\right)}{d p}=-\alpha e^{\alpha h^{2}}\int_{0}^{\infty}e^{-p\alpha\left(u^{2}+h^{2}\right)}\cos\left(ux\right)du=-\frac{1}{2}\sqrt{\pi\alpha}\:e^{\alpha h^{2}}\left(\frac{e^{-p\alpha h^{2}-\frac{x^{2}}{4p\alpha}}}{\sqrt{p}}\right)\tag{4}$$This expression must in turn be integrated with respect to $p$ to find $I_{\mathrm{M}}$:$$I_{\mathrm{M}}=I_{\mathrm{M}}\left(1\right)-I_{\mathrm{M}}\left(\infty\right)=-\int_{1}^{\infty}\frac{dI_{\mathrm{M}}\left(p\right)}{dp}dp=\frac{1}{2}\sqrt{\pi\alpha}\:e^{\alpha h^{2}}\int_{1}^{\infty}\frac{e^{-p\alpha h^{2}-\frac{x^{2}}{4p\alpha}}}{\sqrt{p}}dp\tag{5}$$Mathematica should be able to do this integration as well, but apparently chokes because the coefficient of $p$ in the exponential under the integral is not $-1$. We fix this by rescaling $p$ according to $p=q/(\alpha h^2)$, yielding an integral that Mathematica can handle:$$I_{\mathrm{M}}=\frac{\sqrt{\pi}}{2h}e^{\alpha h^{2}}\int_{\alpha h^{2}}^{\infty}\frac{e^{-q-\frac{h^{2}x^{2}}{4q}}}{\sqrt{q}}dq=\frac{\pi}{4h}e^{\alpha h^{2}}\left(e^{hx}\mathrm{erfc}\left(h\sqrt{\alpha}+\frac{x}{2\sqrt{\alpha}}\right)+e^{-hx}\mathrm{erfc}\left(h\sqrt{\alpha}-\frac{x}{2\sqrt{\alpha}}\right)\right)\tag{6}$$Inserting this into (1) and putting $\alpha=Dt$, we arrive finally at:$$I_{\mathrm{Talon44}}=\sqrt{\frac{D}{\pi t}}\,e^{-\frac{x^{2}}{4Dt}}-Dh\,e^{hx+Dth^{2}}\mathrm{erfc}\left(h\sqrt{Dt}+\frac{x}{2\sqrt{Dt}}\right)\tag{7}$$ which proves the result.

 

[Talon44]

@pasmith

I apologize, I misunderstood what you were trying to accomplish with the simpler example – I did not understand why you were asking me to do other inverse transforms, but now it makes sense. In any case, thank you for taking the time to show that. I did not go through the exercise of the integration by parts yet, but I will do it when I get a chance. Still, it is not clear to me how that substitution helps when I also have a sin(ux) or cos (ux) factor in the integrand. For what it's worth, I tried a similar substitution as one of my first attempts at solving the more difficult problem, and I just got a mess in the sine function that I couldn’t figure out how to resolve.

@renormalize

That’s amazing, I’ve never seen anything like it. Tbh I’m not really great with Mathematica, just the basics really. I considered that there might be other ways to get the program to perform my integration, but I would never have come up with that. I will have to study what you did and see if I can adapt it to some of the other problems I need to solve (which to not appear in my table of transforms). Very much appreciated.

 

[anuttarasammyak]

Jordan’s lemma assures it.

My bad. $e^{-u^2 t}$ term diverges around imaginary axis on R contour.

 

[Talon44]

My bad. $e^{-u^2 t}$ term diverges around imaginary axis on R contour.

So then, how do you evaluate the integral along the semicircular contour as R goes to infinity?

Nvm, I see that an alternative shaped (rectangular) contour may be more effective for integrands containing $e^{-x^2}$. E.g., as shown here. I will try it and see.

 

[anuttarasammyak]

Actually, what I'm really trying to do is show that

2hDπ∫0∞e−Du2tusin⁡uxu2+h2du+2Dπ∫0∞e−Du2tu2cos⁡uxu2+h2du

Has the analytical solution

Dπte−x2/4Dt−hDehx+Dth2erfc(x2Dt+hDt)

The latter expression comes from a table of Laplace transforms.

It may help you to rewrite your integral with non dimensional paratemers, i.e.
$$\sqrt{\frac{D}{t}}\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{e^{-U^2}(HU\sin UX+U^2 \cos UX)}{U^2+H^2}dU$$
where non dimensional parameters are
$$H\equiv h\sqrt{Dt}$$
$$U\equiv u\sqrt{Dt}$$
$$X \equiv \frac{x}{\sqrt{Dt}}$$
Further by suggestion of @renormalize #18
$$-\sqrt{\frac{D}{t}}\frac{1}{\pi}(H\frac{\partial}{\partial X}+\frac{\partial^2}{\partial X^2})\int_{-\infty}^{+\infty}\frac{e^{-U^2}cos UX}{U^2+H^2}dU$$
$$=-\sqrt{\frac{D}{t}}\frac{1}{\pi}(H\frac{\partial}{\partial X}+\frac{\partial^2}{\partial X^2})[e^{-(\frac{X}{2})^2}Re\int_{-\infty}^{+\infty}\frac{e^{-(U\pm \frac{X}{2}i)^2}}{U^2+H^2}dU]$$
It may help you to find an appropriate integral contour.
In an interpretation that D is diffusion constant and t is time, this intefral has dimension of velocity or speed.
The expected transformation is
$$\sqrt{\frac{D}{t}}\ [\ \frac{e^{-(\frac{X}{2})^2}}{\sqrt{\pi}}-He^{HX+H^2}erfc(\frac{X}{2}+H)\ ]$$
$$=\sqrt{\frac{D}{t}}\ \frac{1}{\sqrt{\pi}}\ e^{-(\frac{X}{2})^2}(\ 1-2He^{(\frac{X}{2}+H)^2}\int_0^{\frac{X}{2}+H} e^{-s^2}ds\ )$$
for x,H >0

My attempt, say X,H>0, is
$$\int_{-\infty}^{+\infty}\frac{e^{-(U-\frac{X}{2}i)^2}}{U^2+H^2}dU=\frac{\pi}{H}e^{(H-\frac{X}{2})^2}H_{1/2}(\frac{X}{2}-H)+\int_{-\infty}^{+\infty}\frac{e^{-s^2}}{s^2+H^2-(\frac{X}{2})^2+iXs}ds$$

where H_1/2 is Heviside step function.
$$s^2+H^2-(\frac{X}{2})^2+iXs=(s-Hi+\frac{X}{2}i)(s+Hi+\frac{X}{2}i)$$. I hope this makes sense.