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.

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, 
  • $\begingroup$ If I change the Method to "FiniteElement", it complains about your first boundary condition and states the NeumannValue should be used. If I simply comment out the first condition out, and use Table to display values for the solution for x,t values over the range, I get 1 everywhere. So it appears there's a problem with the first bc. Also, if I change the Method to "Adams" or "BDF" etc., I get messages indicating Differential-Algebraic equations aren't supported by these methods. So NDSolve isn't happy with that first Derivative[1,0][u][0,t] == 0 equation. (MethodOfLines is default BTW.) $\endgroup$
    – user87932
    Nov 22, 2023 at 20:42

3 Answers 3


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

enter image description here

Aha, seems that it's indeed the case.


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:


enter image description here

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


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.

  • $\begingroup$ Well, the explanation for "the reason for NDSolve claiming that you're dealing with a DAE system" isn't quite accurate. (I won't say it's incorrect, though. It's indeed related to the structure you've noticed. ) The essential reason is that after discretizing the D[u[x,t],x, t] term via TensorProductGrid method, current implementation of MethodOfLines isn't clever enough to transform the system to the required standard form for ODE solver. I've talked about the topic a bit here: mathematica.stackexchange.com/a/184285/1871 $\endgroup$
    – xzczd
    Dec 3, 2023 at 0:27
  • $\begingroup$ "while u[x,t] = 1 is a possible solution, is it the only solution, or just a specific one out of a larger family?" It's the only solution, of course. You've forgotten to take the i.c. and b.c.s into consideration. $\endgroup$
    – xzczd
    Dec 3, 2023 at 0:46
  • $\begingroup$ Careful when using Method->FEM for time dependent problems. This will use s spatial discretization for the time variable. $\endgroup$
    – user21
    Dec 3, 2023 at 2:07

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?

  • 1
    $\begingroup$ Hi Roland, thanks for your message, but as I mentioned in the original question, the solution of the problem is already known to be u(x,t)=1. The question is specifically oriented towards the behavior of NDSolve; namely, why it fails to recognize that the solution is u(x,t) = 1. Let me know if this clarifies. $\endgroup$
    – phonon
    Nov 22, 2023 at 17:31
  • $\begingroup$ Just to add a little more detail, running my code you will find that the interpolating solution is stable for about 0 < t < 0.1 and then starts to blow up. Presumably this is not considered to be reasonable behavior? $\endgroup$
    – phonon
    Nov 22, 2023 at 17:34
  • $\begingroup$ The equation is none of elliptic, hyperbolic or parabolic type and so its not possible to ascribe boundary condtions in such a way that make it a "well defined problem" in the sense of Cauchy. This fact correxponds to the existence and uniqueness theorems and the convergency of a lattice approximation for PDEs of the second order. $\endgroup$
    – Roland F
    Nov 22, 2023 at 19:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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