# Trying to solve transient semi-infinite 1D Fick's Law with decay, getting constant for answer

I'm trying to find the solution to $$\frac{\partial f(x,t)}{\partial t} = a\frac{\partial^2 f(x,t)}{\partial x^2} - bf(x,t)$$ where $$f(0,t)=1, \lim_{x\to\infty}f(x,t)=0,f(x,0)=0.$$ I tried:

pde = D[f[x, t], t] == a D[f[x, t], {x, 2}] - b f[x, t];
DSolveValue[{pde, f[0, t] == 1, f[x, 0] == 0, f[Infinity, t] == 0},
f[x, t], {x, t}, Assumptions -> {a > 0, b > 0}];


Mathematica just returns 1, with no errors or warnings. I'm not sure how to make sense of this answer, or what I'm doing wrong. Without the decay term $$bf(x,t)$$, Mathematica gives me the expected solution $$f(x,t)=\text{erfc}(\frac{x}{2\sqrt{at}})$$.

I realize this is probably a well-known problem, but I'm turning to Mathematica because I don't know how to solve it myself or how to effectively search for the solution. Independently of how to do this with Mathematica, I'd be grateful for any tips on what the solution actually is.

• The conditions: f[0, t] == 1, f[x, 0] == 0 contradict each other: what is the value at f[0,0]? Oct 23 '20 at 18:50
• Thank you. I replaced f[x, 0] == 0 with f[x, 0] == If[x != 0, 0, 1] to fix this. But Mathematica gave the same results: The expected Erfc function without decay, and "1" with decay. Oct 23 '20 at 21:48
• Thank you. NDSolve works fine, but I'm looking for an analytical solution. An example physical interpretation would be: On one side of a plane, maintain a fixed concentration of radon, then allow the radon to diffuse into the initially empty space on the other side. You say smooth i.c. and Neumann b.c. are needed, but DSolveValue works without either in the no-decay case. In fact, discontinuous i.c. and Dirichlet b.c. are used in Wolfram's heat flow example in the DSolveValue documentation. Oct 24 '20 at 17:54
• @DanielHuber The discontinuity isn't a problem, and is actually typical of heat-transfer problems where temperatures are raised instantaneously. As for the second derivative in $x$, we have two boundary conditions, one at $x = 0$ and one at infinity. Nov 4 '20 at 7:31
• @yawnoc Hi, it really looks like a bug to me. Even with Assumptions -> {a > 0, b == 0} you get the wrong answer. Please report this to support@wolfram,com. Nov 4 '20 at 9:08

The problem can be solved analytically with Fourier sine transform.

Notice Fourier sine transform has the following property:

$$\mathcal{F}_t^{(s)}\left[f''(t)\right](\omega)=-\omega^2 \mathcal{F}_t^{(s)}[f(t)](\omega)+\sqrt{\frac{2}{\pi }} \omega f(0)$$

as long as $$f(\infty)=0$$ and $$f'(\infty)=0$$. So we apply it to the PDE in $$x$$ direction:

fst[(h : List | Plus | Equal)[a__], t_, w_] := fst[#, t, w] & /@ h[a]
fst[a_ b_, t_, w_] /; FreeQ[b, t] := b fst[a, t, w]
fst[a_, t_, w_] := FourierSinTransform[a, t, w]

neweq = fst[{f[x, 0] == 0, D[f[x, t], t] == a D[f[x, t], {x, 2}] - b f[x, t]}, x, w] /.
f[0, t] -> 1 /. HoldPattern@FourierSinTransform[a_, __] :> a /. f -> (F[#2] &)
(*
{F[0] == 0, F'[t] == -b F[t] - a w (-Sqrt[(2/π)] + w F[t])}
*)

tsol = F[t] /. First@DSolve[neweq, F, t]
(* (a E^(-b t - a t w^2) (-1 + E^(t (b + a w^2))) Sqrt[2/π] w)/(b + a w^2) *)


The last step is to transform back, but sadly InverseFourierSinTransform cannot handle tsol. Still, we've obtained an analytic solution involving integral: $$f(x,t)=\sqrt{\frac{2}{\pi }}\int_0^{\infty } \frac{ a \omega \left(1-e^{-t \left(a \omega ^2+b\right)}\right)}{a \omega ^2+b} \sin (x \omega ) \, d\omega$$

dat = ParallelTable[{x, t,
Sqrt[2/Pi] NIntegrate[tsol Sin[w x] /. a -> 1 /. b -> 1, {w, 0, Infinity},
Method -> "LevinRule", MaxRecursion -> 30] // Quiet}, {t, 0.001, 1, 0.05}, {x,
0.001, 2, 0.1}]; // AbsoluteTiming

Flatten[dat, 1] // ListPlot3D


The result is consistent with yawnoc's series solution:

func = Interpolation@Flatten[dat, 1]

Manipulate[Plot[{func[x, t], fSol[10][x, t]}, {x, 0.001, 1.9}, PlotRange -> {0, 1},
PlotStyle -> {Automatic, {Red, Dashed}}], {t, 0.001, 0.9}]


• +1 This is a much more systematic approach than my ad hocery (in fact I have never properly learned how to do problems where there is more than one similarity variable). Is it possible to show that our answers are indeed the same? Nov 4 '20 at 7:51
• @yawnoc Of course, see the new added GIF. (There exists tiny discrepancy when t is close to 0, it's mainly interpolation error. ) Nov 4 '20 at 8:18
• @yawnoc A strict proof showing the 2 solutions are mathematically equivalent is much harder and beyond my reach, though… Nov 4 '20 at 8:31
• Thank you for this answer! I'm not good with Fourier transforms, but this gives me something to study. And it's good to see the integral form of the answer because it's easier to compare it to the error function solution for the no-decay case. Nov 6 '20 at 18:37
• @coriander As to integral transform, you may refer to Integral Transforms and Their Applications by Lokenath Debnath and Dambaru Bhatta. Nov 6 '20 at 18:45

When I run OP's code on 12.0.0 for Linux x86 (64-bit) (April 7, 2019) it doesn't evaluate, and if I delete the decay term it returns the erroneous 1 - {sum of identically zero terms}. It appears Mathematica tries to find a separation solution in $$x$$ and $$t$$, which doesn't work for a similarity problem such as this.

Indeed we have a similarity problem with dimensionless similarity variables $$\xi = x / \sqrt{a t}$$ and $$\eta = b t$$. Since the dependent variable $$f$$ is dimensionless, we must have $$f = f(\xi, \eta)$$.

We seek separated solutions of the form $$U(\xi) V(\eta)$$:

xiVar = x / Sqrt[a t];
etaVar = b t;
fun = u[xiVar] v[etaVar];
eq =
FullSimplify[
D[fun, t] == a D[fun, {x, 2}] - b fun
/. {
x -> x * xi / xiVar,
b -> b * eta / etaVar
}
// Map[# / (u[xi] v[eta]) &]
, t > 0
]


This gives $$2 \eta \cdot \frac{V + V'}{V} = \frac{\xi U' + 2 U''}{U} = k$$ where $$k$$ is the separation constant. Requiring $$U(\xi = \infty) = 0$$ takes care of both the boundary condition at infinity $$f(\infty, t) = 0$$, and the initial condition $$f(x, 0) = 0$$:

removeC = # /. C[_] -> 1 &;
lim = Limit[#, inf -> Infinity] &;
vSol = DSolveValue[eq[[1]] == k, v[eta], eta] // FullSimplify // removeC
uSol = DSolveValue[{eq[[2]] == k, u[inf] == 0}, u[xi], xi] // lim // removeC


\begin{align} V(\eta) &= \mathrm{e}^{-\eta} \eta^{k/2} \\ U(\xi) &= \mathrm{e}^{-\xi^2/4} H_{-1-k}(\xi/2) \end{align} where $$H$$ is HermiteH.

We therefore have $$f = \sum_{k=0}^\infty A_k \mathrm{e}^{-\eta} \eta^{k/2} \mathrm{e}^{-\xi^2/4} H_{-1-k}(\xi/2).$$ (At first I tried an integral, but that resulted in an inverse Laplace transform with Dirac deltas everywhere, effectively a discrete sum.)

To satisfy the boundary condition at $$x = 0$$, we need \begin{align} 1 = f |_{\xi = 0} &= \sum_{k=0}^\infty A_k \mathrm{e}^{-\eta} \eta^{k/2} H_{-1-k}(0) \\ &= \sum_{k=0}^\infty B_k \mathrm{e}^{-\eta} \eta^{k/2} \end{align} where $$B_k = H_{-1-k}(0) A_k$$. Rearranging, we have \begin{align} \mathrm{e}^{\eta} &= \sum_{k=0}^\infty B_k \eta^{k/2} \\ \sum_{i=0}^\infty \frac{\eta^i}{i!} &= \sum_{i=0}^\infty B_{2i} \eta^{i} + \sum_{i=0}^\infty B_{2i + 1} \eta^{i + 1/2} \end{align} and we see that $$B_{2i} = 1 / i!$$ while $$B_{2i + 1} = 0$$.

Therefore:

1 / i! / HermiteH[-1-2i, 0] // FullSimplify


$$A_{2i} = \frac{B_{2i}}{H_{-1-2i}(0)} = \frac{2 ^ {1 + 2i}}{\sqrt{\pi}}$$ with $$A_{2i+1} = 0$$, and we finally have \begin{align} f &= \sum_{i=0}^\infty \frac{2 ^ {1 + 2i}}{\sqrt{\pi}} \mathrm{e}^{-\eta} \eta^i \mathrm{e}^{-\xi^2/4} H_{-1-2i}(\xi/2) \\ f &= \sum_{i=0}^\infty \frac{2 ^ {1 + 2i}}{\sqrt{\pi}} \mathrm{e}^{-b t} (b t)^i \mathrm{e}^{-x^2/(4at)} H_{-1-2i}\left( \tfrac{x}{2 \sqrt{at}} \right). \end{align}

fSol[numTerms_][x_, t_] :=
Sum[
2 ^ (1 + 2i) / Sqrt[Pi]
Exp[-t] t^i
Exp[-x^2 / (4t)] HermiteH[-1-2i, x / (2 Sqrt[t])]
, {i, 0, numTerms}
];

ital[s_String] := Style[s, Italic];
tValues = Range[0, 2, 0.5] // Rest // Prepend[0.1];
Plot[
Table[fSol[10][x, t], {t, tValues}]
// Append @ Exp[-x] (* steady-state solution *)
// Evaluate
, {x, 0, 4}
, AxesLabel -> {ital["x"] Sqrt @ Row[ital /@ {"b", "a"}, "/"], ital["f"]}
, LabelStyle -> Directive[Black, 16]
, PlotLegends -> LineLegend[
tValues // Append[Infinity]
, LegendLabel -> ital["b"] ital["t"]
]
, PlotStyle -> (ConstantArray[Automatic, Length[tValues]] // Append[Dashed])
]


• Thank you! This is all getting beyond my math skills (hence my turning to Mathematica, of course). So I accepted xzczd's answer based on your own praise for it, but it is helpful for me to see both solutions and to have some idea of the route to get to them. (And just knowing the phrase "similarity problem" is very helpful.) Nov 6 '20 at 18:27
• You're welcome! This is my new favourite question by the way; the above 'method' is in fact me winging it; I've never seen a problem done like this actually. I've also been told that this problem isn't technically a similarity problem, because there's a length scale $\sqrt{a/b}$ and a time scale $1/b$. But the same ideas worked. For an actual similarity problem, see this tutorial of mine. Nov 7 '20 at 3:33