Incorrect Result from DSolve for Diffusion IBVP?

Question

The result of DSolve for the following IBVP doesn't seem to satisfy the IBVP:

\begin{cases} u_{t} - ku_{xx} = 0 &\mbox{} k>0,x>0,t>0 \\ u(0,t)=p(t) &\mbox{} t>0\\ u(x,0)=0 &\mbox{} x>0 \end{cases}

eq = D[u[x, t], t] - k*D[u[x, t], x, x] == 0
ic = u[x, 0] == 0
bc = u[0, t] == p[t]
FullSimplify[DSolve[{eq, ic, bc}, u[x, t], {x, t}, Assumptions -> k > 0 && x > 0 && t > 0]]

This is what I got:

Answer 1

Solution given by DSolve is correct, it just can't be verified by naive substitution.

This problem is similar to, but a bit more involved than your previous one. First of all, as done in my previous answer, we introduce a positive $\epsilon$ to the solution:

eq = D[u[x, t], t] - k D[u[x, t], x, x] == 0;
ic = u[x, 0] == 0;
bc = u[0, t] == p[t];
sol = 
  u[x, t] /. 
    First@DSolve[{eq, ic, bc}, u[x, t], {x, t}, 
        Assumptions -> k > 0 && x > 0 && t > 0]

solfuncmid[x_, t_] = 
 Inactivate[
  sol /. h_[a__, Assumptions -> _] :> h[a] /. {K[2], 0, t} -> {K[2], 
      0, t - ϵ} // Evaluate, Integrate]

Remark

The rule h_[a__, Assumptions -> _] :> h[a] removes the Assumptions option to make the solution look good and avoid unnecessary trouble in subsequent verification, the Inactivate[…] is necessary for v12.0.1 to make subsequent calculation faster, because the Integrate in output of DSolve in v12.0.1 isn't wrapped by Inactive.

Substitute it back to the PDE and combine the integrals:

residual = eq[[1]] /. u -> solfuncmid // Simplify

residual2 = With[{int = Inactive@Integrate}, residual //. 
        HoldPattern[coef1_. int[expr1_, rest_] + coef2_. int[expr2_, rest_]] :> 
          int[coef1 expr1 + coef2 expr2, rest]] // Simplify // Activate

Remark

The . in coef1_. is the shorthand for Optional, it's added so the following type of pattern matching will happen:

aaa /. coef_. aaa -> (coef + 1) b 
(* 2 b *)

Just the same as in the previous answer, when $\epsilon \to 0$ the … Exp[-(…)^2] can be replaced with a … DiracDelta[…]:

residual3 = residual2 /. Exp[coef_ a_^2] :> DiracDelta[a]/Sqrt[-coef] Sqrt[Pi]
(* (x DiracDelta[x] p[t - ϵ])/(Sqrt[k] Sqrt[1/(k ϵ)] ϵ^(3/2)) *)

Given that $x>0$, DiracDelta[x] == 0, so we've verified the solution satisfies the PDE.

Remark

Though Simplify can be used in last step to show residual3 == 0, I've avoided it because of the issue mentioned here.

Verification of initial condition (i.c.) is trivial:

solfuncmid[x, t] /. {t -> 0, ϵ -> 0} // Activate

What's really new compared with the previous problem is the verification of boundary condition (b.c.). The solution only satisfies the b.c. when $x \to 0^+$, so a direct substitution won't work, and doesn't make sense actually, because generally the integral in sol diverges at $x=0$. (Notice Integrate[1/(t - s)^(3/2), {s, 0, t}] diverges. )

Remark

One can also turn to numeric calculation to convince oneself. Here's a quick test with $p(t)=t$:

With[{int = Inactive[Integrate]}, 
  solfuncmid[x, t] /. 
    coef_ int[a__] :> int[a] /. {k -> 1, ϵ -> 0, t -> 2, 
    Integrate -> NIntegrate, x -> 0, p -> Identity}] // Activate
(* NIntegrate::ncvb *)    
(* 2.6163*10^33 *)

To verify the b.c., we transform the solution based on integration by parts:

soltransformed = 
 With[{int = Inactive[Integrate]}, 
  Assuming[{t > K[2], k > 0, x > 0, t > 0, ϵ > 0}, 
    solfuncmid[x, t] /. 
      int[expr_ p[v_], rest_] :> 
       With[{i = Integrate[expr, K[2]]}, 
        Subtract @@ (i p[K[2]] /. {{K[2] -> t - ϵ}, {K[2] -> 0}}) - 
         int[i p'[K[2]], rest]] // Simplify] //.
   coef_ int[a_, b__] :> int[coef a, b]]

Then we take the limit $\epsilon \to 0^+$. It's a pity Limit can't handle soltransformed all at once (this is reasonable of course, the unknown function p[t] is on the way), but by calculating

Limit[Gamma[1/2, x^2/(4 k ϵ)], ϵ -> 0, 
 Direction -> "FromAbove", Assumptions -> {k > 0, x > 0}]
(* 0 *)

separately, we know the correct limit (assuming p[t] is nice enough) is

sollimit = soltransformed /. x^2/(4 k ϵ) -> Infinity /. ϵ -> 0

Now we can substitute $x=0$:

sollimit /. x -> 0 // Activate // Simplify

Integrate refuses the calculate further, which is again reasonable, but it's clear the expression above simplifies to p[t] assuming p[t] is a nice enough function, so the b.c. is verified.

Tested on v12.0.1, v12.1.0.

---

Just for fun, here's a solution based on Fourier sine transform:

Clear@fst
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]

tset = fst[{eq, ic}, x, w] /. Rule @@ bc /. 
  HoldPattern@FourierSinTransform[a_, __] :> a

tsol = DSolve[tset, u[x, t], t][[1, 1, -1]]

The last step is to transform back. Assuming $p(t)$ is a nice enough function so that the order of integration can be interchanged:

With[{int = Inactive[Integrate]}, 
 solfourier = tsol /. 
   coef_ int[a_, rest_] :> 
    int[InverseFourierSinTransform[coef a, w, x], rest]]

It's clear solfourier is equivalent to sol given that $k>0$. Solution verified, once again.

Answer 2

I have not looked too carefully at your simplifications, since hard to read, and it is better to say explicitly why you think the solution is wrong, instead of just showing code, since I was not sure what you are doing there.

The normal way to verify solution from DSolve is to do pde=.... then sol=DSolve[...,u,.....] then pde/.sol//Simplify but this does not simplify to True here, since it does not know what to do with the integral inside.

But I verified by hand that Mathematica's solution is correct.

This is what Mathematica gives as solution

Clear["Global`*"]
eq = D[u[x, t], t] - k*D[u[x, t], x, x] == 0
ic = u[x, 0] == 0
bc = u[0, t] == p[t]
sol = First@DSolve[{eq, ic, bc}, u[x, t], {x, t}, 
    Assumptions -> k > 0 && x > 0 && t > 0];
(sol = sol /. K[2] -> \[Tau])

In Latex, the above is

$$ \Large u(x,t)\to \frac{x \int _0^t\frac{e^{-\frac{x^2}{4 k t-4 k \tau }} p(\tau )}{(t-\tau )^{3/2}}d\tau }{2 \sqrt{\pi } \sqrt{k}} $$

---

And this is my hand solution which gives same answer.

$$ \begin{cases} u_{t} = ku_{xx} &\mbox{} k>0,x>0,t>0 \\ u(0,t)=p(t) &\mbox{} t>0\\ u(x,0)=0 &\mbox{} x>0 \end{cases} $$

Let $U\left( x,s\right) $ be the Laplace transform of $u\left( x,t\right) $ defined as $ \mathcal{L} \left( u,t\right) =\int_{0}^{\infty}e^{-st}u\left( x,t\right) dt$. Applying Laplace transform to the above PDE gives

$$ sU\left( x,s\right) -u\left( x,0\right) =kU_{xx}\left( x,s\right) $$

But $u\left( x,0\right) =0$, the above simplifies to \begin{align*} sU & =kU_{xx}\\ U_{xx}-\frac{s}{k}U & =0 \end{align*}

The solution to this differential equation is

$$ U\left( x,s\right) =c_{1}e^{\sqrt{\frac{s}{k}}x}+c_{2}e^{-\sqrt{\frac{s}{k} }x} $$

Assuming solution $u\left( x,t\right) $ bounded as $x\rightarrow\infty$ and since $k>0$, then $c_{1}=0$. Hence

\begin{equation} U\left( x,s\right) =c_{2}e^{-\sqrt{\frac{s}{k}}x}\tag{2} \end{equation}

At $x=0,u\left( 0,t\right) =p\left( t\right) $. Therefore $U\left( 0,s\right) = \mathcal{L} \left( p\left( t\right) \right) =P\left( s\right) $. At $x=0$, the above gives

$$ P\left( s\right) =c_{2} $$

Hence (2) becomes

\begin{equation} U\left( x,s\right) =P\left( s\right) e^{-\sqrt{\frac{s}{k}}x}\tag{3} \end{equation}

By convolution, the above becomes

\begin{equation} u\left( x,t\right) =p\left( t\right) \circledast G\left( x,t\right) \tag{4} \end{equation}

Where $G\left( x,t\right) $ is the inverse Laplace transform of $e^{-\sqrt{\frac{s}{k}}x}$ which is $\frac{xe^{\frac{-x^{2}}{4kt}}} {2\sqrt{k\pi}t^{\frac{3}{2}}}$, hence (4) becomes

\begin{align*} u\left( x,t\right) & =p\left( t\right) \circledast\frac{xe^{\frac {-x^{2}}{4kt}}}{2\sqrt{k\pi}t^{\frac{3}{2}}}\\ & =\Large \frac{x}{2\sqrt{k\pi}}\int_{0}^{t}\frac{p\left( \tau\right) }{\left( t-\tau\right) ^{\frac{3}{2}}}e^{\frac{-x^{2}}{4k\left( t-\tau\right) } }d\tau \end{align*}

Which as you can see the same as Mathematica solution.

Answer 3

The Mathematica documentation for 12 onward shows an application for the built-in Convolve:

"Obtain a particular solution for a linear ordinary differential equation using convolution:"

First@Simplify[y[x] /. DSolve[{y'[x] + y[x] == x^2, y[0] == 0}, y, x]]

2 - 2 E^-x - 2 x + x^2


Expand[Simplify[Convolve[E^(-t) UnitStep[t], t^2 UnitStep[t], t, x], 
  x > 0]]

2 - 2 E^-x - 2 x + x^2

So this is the new methodology instead of using DSolve to solve ordinary linear differential equations.

In the presentation for Mathematica 11 there was and is still accessible an example solving the heat equation by convolution:solve the boundary value problem of the heat equation.

Hope that helps and clears the adversarial value intention by Wolfram Inc. in the Wolfram Language for the importance of the powerful methodologies of Convolve. Convolve is more elegant and ab-initio than DSolve. It is based on distribution theory and therefore methodologically more advanced modern and numerically stable. It incorporates easier the symbolical solving capabilities of Mathematica and makes Mathematica more superior to competitors at present.

The only drawback is that the example uses a certain initial distribution.

In addition since some former versions of Mathematica this can be solved with GreenFunction. Under the section Examples - basic examples:

"Green's function for the heat operator on the real line:"

GreenFunction[D[u[x, t], t] - D[u[x, t], {x, 2}], 
 u[x, t], {x, -\[Infinity], \[Infinity]}, t, {y, s}]

Plot3D[Evaluate[
  Table[% /. {y -> 1.5, s -> i}, {i, 0.105, 0.5, 0.1}]], {x, 0, 
  3}, {t, 0.1, 0.5}, PlotRange -> {0, 4}]

There is another section in the documentation under Scope concerning Heat Equation.

Green's function for the heat operator in the plane:

GreenFunction[\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[x, y, t]\)\) - Laplacian[
  u[x, y, t],{x, y}], 
 u[x, y, t], {x, y} \[Element] FullRegion[2], t, {m, n, p}]

   GreenFunction[\!\(
    \*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[x, y, t]\)\) - Laplacian[
      u[x, y, t],{x, y}], 
     u[x, y, t], {x, y} \[Element] FullRegion[2], t, {m, n, p}];
    Plot3D[Table[% /. {m -> 0, n -> 0, p -> 0, a -> 1, t -> r}, {r, {0.1, 
         0.2, .3}}] // Evaluate, {x, -2, 2}, {y, -2, 2}, PlotRange ->All]

Now what is the main advantage of the Green's function formalism to the DSolve solution? In the daily usual business none. Both give in all cases the very same results. But the foundational formalism is rather more modern with Green's functions and more general it unifies the very many special cases of DSolve and is more standardized cause by a more coherent development. It is somehow a generalization of DSolve methodologies but that is only for experts important.

Since Green's function are more thorough studies on universities it is up to yourself to keep you compatible in discussions. Mathematica itself offers both and makes the Green's function closer relation to modelling by for example the built-in TransferFunctionModel. Everybody who is interesting in modelling complex technical modern system should be able to operate and design with the most modern Mathematica modelling functions and features.

The steps up to the usage of CreateSystemModell have to be mastered to gain advantage of the all-new Mathematica computation power. Do not by shy to tackle full Mathematica capabilities. That is the main problem and possibly error of this answer by DSolve.

Wolfram Inc. is so friendly to pack some knowledge into Mathematica. For this question look at the built-in SystemModelExamples. And there atop start Modellica.