Trying to solve transient semi-infinite 1D Fick's Law with decay, getting constant for answer
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}]
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])
]