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)$ ?


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}) $$

  • 1
    $\begingroup$ This problem seems to hinge on the math, rather than on Mathematica, doesn't it? Wouldn't it be a better fit for math.stackexchange.com? $\endgroup$
    – MarcoB
    Nov 19, 2019 at 14:32
  • $\begingroup$ @MarcoB I guess you are right to some extent. The transform I posted is the final step of a coupled heat transfer problem I have been working on (regrading which I already took some help from MathStack). Getting the inverse Laplace transform of $(1)$ is something which is supposed to be the next step. I was hoping if that could be done using Mathematica. $\endgroup$
    – Avrana
    Nov 19, 2019 at 15:49
  • $\begingroup$ Perhaps this might be considered in terms of the Fourier transform instead of the Laplace transform. Depending on the values of alpha and beta a stationary phase approximation might be appropriate. $\endgroup$
    – mikado
    Nov 19, 2019 at 19:59
  • $\begingroup$ @mikado I am sorry to say, but I was unable to understand the essence of your remark. In any case, I have added an Origins section to my original question to provide some context to the problem. $\endgroup$
    – Avrana
    Nov 20, 2019 at 3:02
  • $\begingroup$ Are you sure you've deduced $(1)$ correctly? We know Laplace transform is a tool for initial value problem, but you seem to use it in $x$ direction where $x \in [0,1]$? $\endgroup$
    – xzczd
    Nov 20, 2019 at 6:30

1 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


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]


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}]

enter image description here

…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 = 
   Block[{α = 1, β = 2, γ = 3, L = 4, d = 5, Tfi = 6}, expr], 

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}] // 


enter image description here

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.


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

  • $\begingroup$ This is extremely well put together. Really appreciate the effort that you had to put in here. I will go through each step separately now and understand them. $\endgroup$
    – Avrana
    Nov 20, 2019 at 12:15
  • $\begingroup$ Can this problem yield to the FFT method ? $\endgroup$
    – Avrana
    Jul 27, 2020 at 2:58
  • 1
    $\begingroup$ @IndrasisMitra No. Just similar your previous question, $\frac{\partial t_c}{\partial y}$ and $\frac{\partial t_h}{\partial x}$ are on the way. Once again, please check the property of finite Fourier transform carefully to understand what type of problem can be solved with it. $\endgroup$
    – xzczd
    Jul 27, 2020 at 3:06
  • $\begingroup$ I apologize for this. I should have been more careful. Probably I got confused. I do have a two-dimensional variant to the said question where I encounter a partio-integral differential equation in one variable $T$ with no $t$ in b.c.(s). $\endgroup$
    – Avrana
    Jul 27, 2020 at 3:41
  • $\begingroup$ I have posted a two-dimensional version of the aforementioned question here which I think should yield to the transform method. I have concluded this looking at the properties you have mentioned in your original FFT post and some similarities to the question you have answered above. $\endgroup$
    – Avrana
    Jul 27, 2020 at 5:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.