# Inverse Laplace transform help- inverse heat conduction

I've been stuck on an inverse Laplace transform for my research. Would be greatly appreciative of any help solving the inverse Laplace transform of

$$\overline{T}=\frac{2}{k_{1}} \frac{1}{s} \left(\alpha_{1} / s\right)^{1 / 2}\left(\frac{1}{(1+\sigma) \exp \left\lbrace a\left(s / \alpha_{1}\right)^{1 / 2}\right\rbrace-(1-\sigma) \exp \left\lbrace-a\left(s / \alpha_{1}\right)^{1 / 2}\right\rbrace}\right)$$

$$\alpha_1$$, $$a$$, $$k$$, $$\sigma$$ are all constants.

To put this into context I've got a thin layer of insulating material on the surface of a thermocouple. Based on a measured temperature I need to determine the surface heat flux. The equation is Eq. 4 of doi: 10.1016/0017-9310(87)90045-7 if this has anyone interested

Edit: Additional information from here on. This solution has been presented in another paper (doi: 10.1115/1.2752188) as $$T_2(t) = \frac{2\left[1+\left( \frac{1-\sigma}{1+\sigma} \right)\right]\sqrt{\alpha_1}}{k_1} \left\lbrace \sqrt{\frac{t}{\pi}} \exp\left({\frac{a^2}{-4\alpha_1 t}}\right) - \left(\frac{a}{2\sqrt{\alpha_1}}\right) \mathrm{erfc}\left(\frac{a}{4\sqrt{\alpha_1} \sqrt{t}}\right) + \sum_{n=1}^{\infty} \left( \frac{1-\sigma}{1+\sigma} \right)^n \left( \sqrt{\frac{t}{\pi}} \exp{\left(\frac{-k_a^2}{4t}\right)} - \frac{k_a}{2} \mathrm{erfc}\left( \frac{k_a}{2\sqrt{t}} \right) \right) \right\rbrace$$ where $$k_a = \frac{(2n+1)a}{\sqrt{\alpha_1}}$$. However, when I use this equation to process my experimental data the match is poor, especially for small $$t$$. If I change to $$k_a = \frac{(2n+0)a}{\sqrt{\alpha_1}}$$ the heat flux profile is perfect.

Because the experimental data matches the $$k_a = \frac{(2n+0)a}{\sqrt{\alpha_1}}$$ model I'm convinced it is correct, not $$k_a = \frac{(2n+1)a}{\sqrt{\alpha_1}}$$ (sample figure link at bottom of post). This is true for tests using square and Gaussian temporal profile heat flux. I've been stuck for over a year on/off on this problem and my thesis submission is coming up, and I'd like to get the work submitted to a journal before then. I'm looking for someone to tell me that my modification is valid, or for someone to call me an idiot if that's the case.

Variables: $$\sigma$$ is positive and greater than 1, in my actual case the value is about 4. $$\alpha_1$$ is positive (about $$118\times 10^{-9}$$), $$a$$ is positive (about $$20\times 10^{-6}$$), and $$0\leq t \leq 1.5$$ seconds.

Sample result of the two models Black- irradiance applied, red- heat flux from an experimental calibration, blue uses $$k_a = \frac{(2n+1)a}{\sqrt{\alpha_1}}$$, green uses $$k_a = \frac{(2n+0)a}{\sqrt{\alpha_1}}$$. The amplitude differences are because the apsorptivity of the material has to be determined which I am doing using results from the experimental (red) and analytical model (green or blue).

Note- more than happy to give credit on any publication that arises (currently with the unmodified model because we can justify using it)

• Do we know if $\sigma$ is smaller or larger than 1? – eyeballfrog May 3 '19 at 23:25
• @Mattos - I believe so, updated post to show what has been presented as a solution before, but it's poor at small $t$. Added my "correction" for what I think the solution is. – sleepy_panda May 3 '19 at 23:52
• Does the model with $2n+1$ fit your data if you replace $\sigma$ with $1/\sigma$? If so I think I know what's going on here. – eyeballfrog May 4 '19 at 0:18
• @eyeballfrog sorry, $\sigma > 1$ (early morning error) . No, does not fit model if changed to $1/\sigma$, but I've just found with some tweaking that the $k_a = 2n+1$ model shape is improved when $\sum_{n=0}^{\infty}$, but the result seems to be out by a factor of $\pi$ – sleepy_panda May 4 '19 at 1:03
• I have not looked into the paper, but just from glancing at $\overline{T}$ with kernel ${\rm e}^{st}$ you can transform the integrand with substitution $v=\sqrt{s}$ into a holomorph one. The new contour is a hyperbola with angle ${\rm e}^{\pm i\pi/4}$ at the asymptotics. Probably the denominator can then be converted using the geometric series and the integrals are elementary (error function). Since the integrand is holomorph the contour doesn't matter and you can evalute at $\lim_{R\rightarrow \infty} R{\rm e}^{\pm i\pi/4}$. Have you tried that? – Diger May 4 '19 at 1:17

To solve this Laplace transform factor $$(1+\sigma)\exp(\sqrt{\frac{s}{\alpha_1}})$$ and expand into a geometric series:

$$\bar{T}(s)=\frac{2\sqrt{\alpha_1}}{k_1(1+\sigma)}s^{-3/2}\sum_{n=0}^{\infty}\Big(\frac{1-\sigma}{1+\sigma}\Big)^n\exp\Big(-(2n+1)a\sqrt{\frac{s}{\alpha_1}}\Big)$$

And Mathematica returns for the ILT:

$$F(s)=s^{-3/2}e^{-a\sqrt{s}}\rightarrow f(t)= 2 e^{-a^2/4 t}\frac{\sqrt{t}}{\sqrt{\pi}} - a~\textrm{erfc}(\frac{a}{2\sqrt{t}})]$$

and substituting $$k_s=\frac{a(2n+1)}{\sqrt{\alpha_1}}$$, we get

$$T(t)=\frac{2\sqrt{\alpha_1}}{k_1(1+\sigma)}\sum_{n=0}^{\infty}\Big(\frac{1-\sigma}{1+\sigma}\Big)^n[2\sqrt{\frac{t}{\pi}}\exp(-\frac{k_s^2}{4t})-k_s~\textrm{erfc}(\frac{k_s}{2\sqrt{t}})]$$

which is exactly the proposed solution, modulo elementary algebra.

EDIT: Calculating $$I(t)=\frac{1}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}e^{-a\sqrt{s}}e^{st}dt$$.

Using the closed contour C that connects the lines $$Re(z)=\sigma$$,the circle in the infinite left half plane and a keyhole contour around the essential singularity at $$s=0$$ we obtain:

$$0=\frac{1}{2\pi i}\oint_C e^{-a\sqrt{z}}e^{zt}dz=\frac{1}{2\pi i}(\int_{0}^{\infty}e^{-ia\sqrt{r}-rt}dr-\int_{0}^{\infty}e^{ia\sqrt{r}-rt}dr+\int\limits_{z=Re^{i\theta}, \frac{\pi}{2}\leq\theta\leq\frac{3\pi}{2}}dze^{-a\sqrt{z}}e^{zt}-\int\limits_{z=\epsilon e^{i\theta}, -\pi\leq\theta\leq\pi }dze^{-a\sqrt{z}}e^{zt}+\int_{\sigma-i\infty}^{\sigma+i\infty}e^{-a\sqrt{s}}e^{st}dt)$$.

Taking the limits $$R\rightarrow\infty, \epsilon\rightarrow 0$$ we obtain that

$$I(t)=\frac{1}{\pi}\int_{0}^{\infty}\sin{(a\sqrt{r})}e^{-rt}dt=-\frac{\partial}{\partial a}\int_{-\infty}^{\infty}\cos(ax)e^{-x^2 t}=\frac{a}{2\sqrt{\pi}t^{3/2}}e^{-a^2/4t}$$

The calculation is similar for other powers of z in the denominator, but is more subtle potentially because of the worsening singularity at z=0. However, in the present case one can obtain the desired formula by integrating twice, once in t and once in a and applying appropriate boundary conditions.

• Actually wolframalpha didn't find a solution, how long did Mathematica calculate? – Diger May 4 '19 at 2:25
• 0.4 secs, for free edition of wolfram alpha the input has to be grossly oversimplified and no free parameters. Do you need a derivation of this Laplace transform? – DinosaurEgg May 4 '19 at 2:27
• I'd be only interested in how to manually calculate $$I(t)=\int \frac{2\sqrt{\alpha_1}}{v^2} \, {\rm e}^{tv^2-\frac{(2n+1)a}{\sqrt{\alpha_1}} \, v} \, {\rm d}v \, .$$ I tried a couple of partial integrations, but that lead nowhere. $I'(t)$ can be simply calculated, but else I dunno. – Diger May 4 '19 at 2:39
• Cheers, you're all legends. Even if you didn't tell me what I wanted to hear! – sleepy_panda May 4 '19 at 2:39
• @Diger You need to use complex analysis or differentiation under the integral sign to compute this integral – DinosaurEgg May 4 '19 at 2:41

An attempt to do this through a slightly different route.

First let's get rid of those constants using the identity $$F(s/\nu)/\nu = \mathcal{L}[f(\nu t)]$$. Then $$T_2(t) = f(\nu t)$$, where $$f$$ satisfies $$\mathcal{L}[f](s) = \frac{2a}{k_1}\frac{1}{s^{3/2}}\frac{1}{(\sigma+1)e^\sqrt{s}+(\sigma-1)e^{-\sqrt{s}}},$$ which now at least only has one parameter to worry about in the transform. Unfortunately, now we have to actually do it. Using the contour $$(\alpha+iu)^2$$ for the inverse Laplace transform and playing around with the resulting integral a bit, I managed to come up with $$f'(t) = \frac{a\sigma}{\pi k_1}\int_{-\infty}^\infty \frac{e^{-tu^2}\cos(u)}{1+(\sigma^2-1)\cos(u)^2}du.$$ Despite its relatively simple form, Mathematica won't do this integral and I didn't find it in Gradshteyn and Rhyzik. Attempting to expand the denominator and make an infinite sum just turns into the infinite sum in the OP, but I'm curious if there is a form for the integral directly.

• I think it can be done in terms of Jacobi theta functions and their integrals/derivatives. – DinosaurEgg May 6 '19 at 21:48