Inverse Laplace transform of this complicated function

Question

I have been solving a coupled PDE system analytically and I need to find the inverse Laplace transform of $(1)$ and get $T(x,y)$. $s$ is the Laplace domain variable and $\alpha, \beta, \gamma, T_{fi}, A , d$ are constants.

$$ \mathcal{L_x}(T)=\frac{As(s+\alpha) - \beta T_{fi}}{s(s^2 - \beta + \alpha s)}+ \frac{\gamma (e^{\sigma y}+e^{\sigma(2d-y)})}{s\sigma (1-e^{2\sigma d})} \tag 1 $$ where

$$\sigma = \sqrt{\frac{\beta s - \alpha s^2 - s^3}{s+\alpha}}$$

I have tried the Wolfram alpha widget but it fails to do the job. The problem comes especially with the second term. Can anyone suggest a way to handle the inverse transformation of $(1)$ ?

---

ORIGINS

I have been trying to solve $(1.1)$, and the Laplace transform mentioned in equation $(1)$ comes from $(1.1)$ subjected to the given boundary conditions.

The term $A$ is $T(0,y)$, which is temporarily being treated as a constant $$ \nabla^2 T - \beta T + \beta\Bigg[\alpha e^{-\alpha x} \Bigg(\int_0^x e^{\alpha s}T(s,y)\mathrm{d}s+\frac{T_{fi}}{\alpha}\Bigg)\Bigg]=0 \tag {1.1} $$ (1.1) is dictated by the following boundary conditions: $$ \frac{\partial T}{\partial x} \vert_{x=0} = \frac{\partial T}{\partial x} \vert_{x=L} = \frac{\partial T}{\partial y} \vert_{y=d} = 0 , \frac{\partial T}{\partial y} \vert_{y=0}=\gamma $$

---

Intermediate steps between $(1.1)$ to $(1)$

Taking Laplace transform of $(1.1)$ w.r.t. $x$ $$ s^2 \mathcal{L_x}T(x,y) - \color{red}{sT(0,y)} - \color{green}{\frac{\partial T(0,y)}{\partial x}} + \mathcal{L_x}\Bigg(\frac{\partial^2 T}{\partial y^2}\Bigg)-\\ \beta \mathcal{L_x}T(x,y) + \frac{\alpha \beta}{\alpha +s}\mathcal{L_x}T(x,y) + \frac{\beta}{\beta +s} T_{fi} = 0 \tag 2 $$ $T(0,y)$ is an unknown and we denote it with the letter $A$ for the rest of this analysis. $$ \mathcal{L_x}\Bigg(\frac{\partial^2 T}{\partial y^2}\Bigg)=\frac{\partial^2}{\partial y^2}\mathcal{L_x}(T(x,y)) $$ Equation $(2)$ becomes $$ \frac{\partial^2}{\partial y^2}\mathcal{L_x}(T)+\Bigg(s^2 - \beta + \frac{\alpha \beta}{\alpha +s}\Bigg)\mathcal{L_x}(T)-sA+\frac{\beta T_{fi}}{\alpha +s}=0 \tag 3 $$ Solve $(3)$ (an O.D.E) to find $\mathcal{L_x}(T)$ $\color{Blue}{\Rightarrow}$ Use $y$ B.C.(s) to determine the constants $\color{Blue}{\Rightarrow}$ Find $\color{black}{T=\mathcal{L_x^{-1}}(T)}$ $\color{Blue}{\Rightarrow}$ Use the $x=L$ B.C. to determine $A$

Solving $(3)$ gives $$ \mathcal{L_x}(T)=\frac{As(s+\alpha) - \beta T_{fi}}{s(s^2 - \beta + \alpha s)}+C_1 e^{\sigma y} + C_2 e^{-\sigma y} \tag 4 $$ where, $$\sigma = \sqrt{\frac{\beta s - \alpha s^2 - s^3}{s+\alpha}}$$ The $y$ boundary conditions become:

$\frac{\partial T}{\partial y} \vert_{y=d} = 0 , \frac{\partial T}{\partial y} \vert_{y=0}=\gamma \color{Blue}{\Rightarrow} \frac{\partial \mathcal{L_x}(T)}{\partial y} \vert_{y=d} = 0,\frac{\partial \mathcal{L_x}(T)}{\partial y} \vert_{y=0}=\gamma$

Utilizing these conditions $C_2=C_1 e^{2\sigma d}$ and $C_1=\frac{\gamma}{s\sigma (1-e^{2\sigma d})}$

Substituting in $(4)$ gives us $$ \mathcal{L_x}(T)=\frac{As(s+\alpha) - \beta T_{fi}}{s(s^2 - \beta + \alpha s)}+ \frac{\gamma (e^{\sigma y}+e^{\sigma(2d-y)})}{s\sigma (1-e^{2\sigma d})} \tag 1 $$

---

For separation of variables I assumed the following ansatz

$$ T(x,y)=\sum_{k=0}^{\infty}f_k(y)\cos(\frac{k\pi x}{L})=f_0(y)+\sum_{k=1}^{\infty}f_k(y)\cos(\frac{k\pi x}{L}) $$

Answer 1

This post contains several code blocks, you can copy them easily with the help of importCode.

---

As already mentioned in the comment above, the deduction of $(1)$ is incorrect because OP forgot $A$ cannot be treated as constant when solving ODE $(3)$, so it doesn't make much sense to continue discussing the Laplace inversion of $(1)$. Since OP's target is just to solve $(1.1)$ analytically, I'll show my solution based on finite Fourier cosine transform and its inversion as an answer. The code is a bit advanced, please check the document carefully by pressing F1 to understand it.

We first interpret the PDE and b.c.s to Mathematica code:

eq = Laplacian[
    T[x, y], {x, y}] - β T[x, 
     y] + β (α Exp[-α x] (Integrate[
         Exp[α s] T[s, y], {s, 0, x}] + Tfi/α)) == 0

bcx = {D[T[x, y], x] == 0 /. x -> 0, D[T[x, y], x] == 0 /. x -> L}

bcy = {D[T[x, y], y] == γ /. y -> 0, D[T[x, y], y] == 0 /. y -> d}

It's easy to notice the integral inside eq can be eliminated:

neweq = eq /. Solve[D[eq, x], Integrate[E^(α s) T[s, y], {s, 0, x}]][[1]] // 
  Simplify[#, α != 0] &

$$\alpha \frac{\partial^2 T}{\partial y^2}+\alpha \frac{\partial^2 T}{\partial x^2}+\frac{\partial^3 T}{\partial y^2 \partial x}+\frac{\partial^3 T}{\partial x^3}=\beta \frac{\partial T}{\partial x}$$

The differential order in $x$ direction becomes $3$, so we need one more b.c., this can be deduced by setting $x$ to $0$ in eq:

newbc = eq /. x -> 0

OK, let's begin solving. Definition of finiteFourierCosTransform and inverseFiniteFourierCosTransform isn't included in this post, please find them in the link above. We make finite Fourier cosine transform in the range $y \in [0, d]$:

rule = finiteFourierCosTransform[a_, __] :> a;

tneweq = finiteFourierCosTransform[neweq, {y, 0, d}, n] /. 
   Rule @@@ Flatten@{bcy, D[bcy, x]} /. rule

tbcx = finiteFourierCosTransform[bcx, {y, 0, d}, n] /. rule

tnewbc = finiteFourierCosTransform[newbc, {y, 0, d}, n] /. (Rule @@@ bcy /. x -> 0) /. rule

Remark

I've stripped off finiteFourierCosTransform because DSolve has difficulty in understanding expression like finiteFourierCosTransform[T[x, y], {y, 0, d}, n]. Just remember that T[x, y] actually denotes finiteFourierCosTransform[T[x, y], {y, 0, d}, n] in tneweq, tbcx and tnewbc.

{tneweq, tbcx, tnewbc} forms a boundary value problem of ODE, it can be easily solved by DSolve:

tsolzero = T[x, y] /. 
  First@DSolve[Simplify[#, n == 0] &@{tneweq, tbcx, tnewbc}, T[x, y], x]

tsolrest = T[x, y] /. 
  First@DSolve[Simplify[#, n > 0] &@{tneweq, tbcx, tnewbc}, T[x, y], x]

tsol = Piecewise[{{tsolzero, n == 0}}, tsolrest]

Remark

The n == 0 case is solved separately, or DSolve won't handle the removable singularity properly.

The final step is to transform back:

sol = inverseFiniteFourierCosTransform[tsol, n, {y, 0, d}]

…As already mentioned, the solution is rather complicated.

"So, how do you know the mess is correct?" OK, let's verify it by solving the problem numerically. However, the somewhat strange newbc stops us from using NDSolve, so I'll solve the problem based on FDM. I'll use pdetoae for the generation of finite difference equations:

setparameters = 
  Function[expr, 
   Block[{α = 1, β = 2, γ = 3, L = 4, d = 5, Tfi = 6}, expr], 
   HoldAll];

test = Compile[{x, y}, #] &[sol /. C -> 20 // ReleaseHold // ToRadicals] // setparameters;

points@x = points@y = 50; domain@x = {0, L}; domain@y = {0, d};
(grid@# = Array[# &, points@#, domain@#]) & /@ {x, y};
difforder = 2;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[T[x, y], grid /@ {x, y}, difforder];

delx = #[[3 ;; -2]] &; dely = #[[2 ;; -2]] &;

ae = dely /@ delx@ptoafunc@neweq;
aebcx = dely /@ ptoafunc@bcx;
aebcnew = dely@ptoafunc@newbc;
aebcy = ptoafunc@bcy;
var = Outer[T, grid@x, grid@y] // setparameters // Flatten;
{barray, marray} = 
 CoefficientArrays[Flatten@{ae, aebcx, aebcnew, aebcy} // setparameters, var]

nsol = ListInterpolation[Partition[LinearSolve[marray, -N@barray], points@y], 
   grid /@ {x, y}] // setparameters

 lst = Table[
    Plot[{test[x, y], nsol[x, y]}, {x, 0, L}, 
     PlotLegends -> {"Series Solution", "FDM Solution"}], {y, 0, d, d/25}] // 
   setparameters;

ListAnimate@lst

As we can see, the 2 solutions agree well, and will be better if you increase points[x], points[y], number of terms in test.

Remark

The discrepancy at $y=0$ and $y=d$ is slightly large, this is expected, because the b.c.s are actually inconsistent.