NDSolve fails to solve trivial PDE with mixed partial derivatives

Question

Consider the following PDE with independent variables $(x,t) \in [0,1]\times[0,\infty)$ $$2u_{xt}(x,t)=u_{xx}(x,t)$$ initial condition $u(x,0)=1$ and with boundary conditions $u_x(0,t)=0$ and $u(1,t) =1$.

Despite the fact that the solution is clearly $u(x,t)=1$, NDSolve seems unable to yield stable results. I've tried increasing WorkingPrecision and changing to method of lines to no avail.

Clear["Global`*"]
T = 1;

pde = {2 Derivative[1, 1][u][x, t] == Derivative[2, 0][u][x, t]};
bcs = {Derivative[1, 0][u][0, t] == 0, u[1, t] == 1};
ic = {u[x, 0] == 1};
system = Join[pde, ic, bcs]
interp = NDSolveValue[system, u, {x, 0, 1}, {t, 0, T}]

Animate[Plot[interp[x, t], {x, 0, 1}, PlotRange -> {0, 1.5}], {t, 0, 
  T}]

Answer 1

Rule of thumb: if NDSolve doesn't work well on solving PDE, and you're sure you haven't made any simple mistake, then usually the problem lies in spatial discretization.

The following paragraph in the tutorial The Numerical Method of Lines seems to be related:

Fourth-order differences typically provide a good balance between truncation (approximation) error and roundoff error for machine precision. However, there are some applications where fourth-order differences produce excessive oscillation (Gibb's phenomena), so second-order differences are better.

So, let's try turning to second order spatial discretization instead of the default fourth order:

mol[n : _Integer | {_Integer ..}, o_ : "Pseudospectral"] := {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
        "MinPoints" -> n, "DifferenceOrder" -> o}}

interp = NDSolveValue[system, u, {x, 0, 1}, {t, 0, T}, Method -> mol[25, 2]];

Plot3D[interp[x, t], {x, 0, 1}, {t, 0, T}, PlotRange -> {0, 2}]

Aha, seems that it's indeed the case.

Answer 2

I tried changing the Method to see what would happen. I first tried "Adams", since it was first in the list, and was informed that this method couldn't handle a Differential-Algebraic system. I tried a few others and got the same message. I eventually tried Method->"FiniteElement" and got the message ref/message/NDSolveValue/fembderiv, which basically states that the FEM requires a NeumannValue instead of Derivative[..] to handle spatial derivatives. Since the default NeumannValue is zero unless otherwise specified, I commented it out of your set and tried running the following code:

pde = {2 Derivative[1, 1][u][x, t] == Derivative[2, 0][u][x, t]};
bcs = {(*Derivative[1, 0][u][0, t] == 0, *)u[1, t] == 1};
ic = {u[x, 0] == 1};
system = Join[pde, ic, bcs];
interp = Flatten@NDSolveValue[system, u, {x, 0, 1}, {t, 0, 1},Method->"FiniteElement"];

This ran without generating any error messages, and I obtained the following 3D plot:

Plot3D[interp[x,t],{x,0,1},{t,0,1},PlotRange->{All,All,{0,1.2}}]

So this agrees with your expectations. But why did you encounter problems with the default method, and what is the reason for NDSolve claiming that you're dealing with a DAE system? Let's look at your pde, which can be written in this form:

2 D[u[x,t],x, t] - D[u[x,t],x,x] == 0 

This can be expressed as

D[f[x,t],x] == 0

where

f[x,t] = D[u[x,t],t] - D[u[x,t],x] 

Since the derivative of f[x,t] wrt x is zero, it's a constant wrt x, i.e.

f[x,t] = g[t]

for some arbitrary function g. This is why NDSolve classifies your equation as a DAE. Now while u[x,t] = 1 is a possible solution, is it the only solution, or just a specific one out of a larger family? I'm not sure, but either way, this seems a bit subtle, and may explain why you're encountering numerical instability.

One other suggestion. NDSolve doesn't generate messages when run in your original configuration, so instead of comparing the output to what you expected to get, it might be better to see if the result satisfies the equations you input. There are two methods for verifying solutions here: FEMDocumentation/tutorial/FiniteElementBestPractice#754253806

The straightforward way is to plug the solution back into the pde and compare the results to zero on the grid. The second, which is better, is to use the method of "manufactured solutions", which avoids introducing additional numerical errors while verifying the solutions.

Answer 3

Your equation is a differential condition of first order for $u_x$ with solution $$\text{DSolveValue}\left[\frac{\partial ^2u(x,t)}{\partial x\, \partial x}=2 \frac{\partial ^2u(x,t)}{\partial x\, \partial t},u,\{t,x\}\right]\text{/.}\, \{c_1\to a,c_2\to b\}$$

    Function[{x,t},a[t]+b[t+2 x]]

leaving us with the equations

    {a[0] + b[2 x] == 1, 2 Derivative[1][b][t] == 0, a[t] + b[2 + t] == 0}

Any ideas?