5
$\begingroup$

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.

$\endgroup$
7
  • 2
    $\begingroup$ The conditions: f[0, t] == 1, f[x, 0] == 0 contradict each other: what is the value at f[0,0]? $\endgroup$ Oct 23 '20 at 18:50
  • $\begingroup$ 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. $\endgroup$
    – coriander
    Oct 23 '20 at 21:48
  • 1
    $\begingroup$ 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. $\endgroup$
    – coriander
    Oct 24 '20 at 17:54
  • 2
    $\begingroup$ @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. $\endgroup$
    – yawnoc
    Nov 4 '20 at 7:31
  • 1
    $\begingroup$ @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. $\endgroup$ Nov 4 '20 at 9:08
6
$\begingroup$

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

enter image description here

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

enter image description here

$\endgroup$
5
  • $\begingroup$ +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? $\endgroup$
    – yawnoc
    Nov 4 '20 at 7:51
  • 1
    $\begingroup$ @yawnoc Of course, see the new added GIF. (There exists tiny discrepancy when t is close to 0, it's mainly interpolation error. ) $\endgroup$
    – xzczd
    Nov 4 '20 at 8:18
  • $\begingroup$ @yawnoc A strict proof showing the 2 solutions are mathematically equivalent is much harder and beyond my reach, though… $\endgroup$
    – xzczd
    Nov 4 '20 at 8:31
  • $\begingroup$ 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. $\endgroup$
    – coriander
    Nov 6 '20 at 18:37
  • 1
    $\begingroup$ @coriander As to integral transform, you may refer to Integral Transforms and Their Applications by Lokenath Debnath and Dambaru Bhatta. $\endgroup$
    – xzczd
    Nov 6 '20 at 18:45
5
$\begingroup$

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

Solution plot

$\endgroup$
2
  • $\begingroup$ 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.) $\endgroup$
    – coriander
    Nov 6 '20 at 18:27
  • 1
    $\begingroup$ 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. $\endgroup$
    – yawnoc
    Nov 7 '20 at 3:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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