# What can one do with extremely stiff problem in NDSolve?

Consider the following illustrative problem:

$$\frac {\partial f} {\partial t} = \frac {\partial} {\partial x}(x f) + \frac {\partial} {\partial x}(f \frac {\partial f} {\partial x})$$

This is some nonlinear advection-diffusion problem, which after some time will come to a steady-state. The objective is to investigate the this steady-state properties, but unfortunately the problem becomes very stiff and the equilibrium is not reached. My solution in Mathematica:

time = .9;  (* solution stop time *)

y0[x_] := Exp[-x^2/2];  (* initial distribution *)
eq = D[y[x, t], t] == D[x y[x, t], x] + D[y[x, t] D[y[x, t], x], x]; (* PDE to solve *)
bc = {y[x, 0] == y0[x], y[-6, t] == y[6, t] == 0};  (* boundary conditions *)

sol = NDSolveValue[{eq, bc}, y, {x, -6, 6}, {t, 0, time}, AccuracyGoal->5, PrecisionGoal->5,
Method -> {"MethodOfLines", Method->"StiffnessSwitching",
"SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->9, "MinPoints"->3000}}];


After some time the distribution has very sharp edges and solving stops:

But the steady-state will take a place only closer to time=3 or so, look at the evolution of the distribution center:

If I increase time even to 1 the solution immediately blows-up. Is there anything one can do in Mathematica to get closer to the final stead-state? Many thanks for any help!

• (1) I would start with an investigation of the steady state solution. In a symmetric situation (f'[0]==0) there are two solutions, f1 = 0, and f2 = - x^2/2 + const. None of them shows a non trivial and finite profile. (2) strictly speaking, your initial and boundary conditions are not compatible. (3) please describe where the specific nonlinar term comes from. Apr 17, 2016 at 15:34
• If you're only interested in steady state, is there a reason why you're not you only solving in space? Apr 17, 2016 at 15:45
• @MathX well, yes. Besides steady-state distribution I need also to know when the system will come to an equilibrium.. so I need to study the evolution. In addition, if you solve only steady state problem you - as boundary you only have the initial area (integral of initial distribution), so this becomes rather integro-differential equation, which can easily be solved with some tricks, but.. And the last one - I often encounter this types of problems in PDE, so it's very interesting to know whether there are any techniques to handle this. Sorry not clarified all this in question. Apr 18, 2016 at 8:01
• @Dr.WolfgangHintze (1) please see my previous comment. (2) you're right, but this is beyond solving accuracy, nevertheless to be on the safe side one can multiply the initial distribution with UnitBox or make a piecewice func. (3) not clear what you mean. By nonlinearity here I understand non-const coefficients under derivatives. In reality there are quite comlicated functions of f and x. Apr 18, 2016 at 8:06
• @funnypony I have started a study of your interesting problem, but I have to aks for some patience to write it down. Apr 18, 2016 at 9:43

EDIT #2

My error was useful. It brought me to the conclusion that the difficulties in solving the PDE of the OP are due to the drift term

$$\frac{\partial (x u(x,t))}{\partial x}$$

If the drift term is included, many boundary problems are ill defined. It turns out that there are cases where mathematically there is only a trivial solution u = 0 but numerically you obtain nice non trivial pictures.

Example (notice the linear dispersion term)

$$\frac{\partial u(x,t)}{\partial t}=\frac{\partial ^2u(x,t)}{\partial x^2}+\frac{\partial (x u(x,t))}{\partial x}$$

$$u(x=\pm 1,t)=0$$

$$u(x,0)=\cos ^2\left(\frac{\pi x}{2}\right)$$

I shall describe the details soon.

EDIT

After some days of consideration I found that with a drastic simplification we can find a satisfactory solution.

The simplification consists in dropping the "drift" term in the PDE.

The "energy leak paradoxon" which persists here, was solved.

Further discussion and revision of the original post are still to be done.

With an intial condition quit close to that of the OP

h[x_] := Exp[-18 x^2]


the solution is (notice the missing "drift" term)

uu[x_, t_] =
u[x, t] /.
NDSolve[D[u[x, t], t] == D[u[x, t] D[u[x, t], x], x] &&
u[-1, t] == 0 && u[+1, t] == 0 && u[x, 0] == h[x],
u[x, t], {x, -1, 1}, {t, 0, 10}][[1]];


and a 3D plot of the function

Plot3D[uu[x, t], {x, -1, 1}, {t, 0, 1}, PlotRange -> All,
AxesLabel -> {"x", "t", "u"}, LabelStyle -> {12, Bold},
ViewPoint -> 10/22 {Pi/2, -3 Pi/2 , \[Pi]},
PlotLabel ->
Style["Numerical solution of PDE\n\
\!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]\)u = \!$$\*SubscriptBox[\(\ \[PartialD]$$, $$x$$]\)(u \!$$\*SubscriptBox[\(\[PartialD]$$, \
$$x$$]\)u)\nu(x=\[PlusMinus]1,t) = 0\nu(x,t = 0) = Exp[-(6 \
x\!$$\*SuperscriptBox[\()$$, $$2$$]\)/2]\n", 14, Bold]]


There is no numeric instability seen here anymore.

This goes also for initial conditions of the form

h1[x_] = Cos[ Pi/2 x ]^k  (* k=1,2,3,4 *)


The stationary state reached for large t is the trivial one, u = 0.

I took some days off in order to clarify the "mystery" of the energy leak which is in fact specific to this non linear problem because the function appears as a factor before the derivative. This is in contrast to the situation with the diffusion equation. The solution - in brief - is that u(x) u'(x) can remain finite for x->1 if the derivative diverges. As an example take f(x) = Sqrt(1-x^2) which gives f(x)f'(x) = - x. The numerical results for large enough t show large (divergent) spatial derivatives at the boundaries. Hence the "energy" will leak out indefinitely until it is zero.

The "energy" of the inital state is

w0 = Integrate[h[x], {x, -1, 1}]

(* Out[13]= 1/3 Sqrt[\[Pi]/2] Erf[3 Sqrt[2]] *)


The normalized "energy" as a function of time is then given by

w[t_] := Integrate[uu[x, t], {x, -1, 1}]/w0


Plots over different regions of time are

Plot[w[t], {t, 0, 1}, PlotRange -> {0, 1.1},
PlotLabel ->
Style["normalized 'energy' as a function of time\ninitial condition u(x,t = \
0) = Exp[-(6 x\!$$\*SuperscriptBox[\()$$, $$2$$]\)/2]" , 14, Bold],
AxesLabel -> {"t", "w"}, LabelStyle -> {12, Bold}]


Original post

The following is the attempt of an in depth study of this interesting problem. Due to a lack of time here are only the first steps stating some of the specifiy problems encountered. I shall complete it subsequently.

The simplified problem

To begin with, let us simplify the problem without losing the main features. We consider the interval -1 < x < 1 instead of the infinite interval, and take as inital condition a power of (1-x^2).

The PDE is

eq = D[f[x, t], t] == D[x f[x, t], x] + D[f[x, t] D[f[x, t], x], x]

(*
Out[6]=
Derivative[0, 1][f][x, t] ==
f[x, t] + x*Derivative[1, 0][f][x, t] + Derivative[1, 0][f][x, t]^2 +
f[x, t]*Derivative[2, 0][f][x, t]
*)


The initial distribution is taken as

f0[x_, k_] = (1 - x^2)^k;


The boundary conditions are

bound = {f[1, t] == 0 && f[-1, t] == 0};


The numerical solution is then (with k = 4, say)

sol = NDSolve[eq && f[x, 0] == f0[x, 4] && bound,
f[x, t], {x, -1, 1}, {t, 0, 20}][[1]];

ff[x_, t_] = f[x, t] /. sol;


The solution approaches the steady state solution in a few time steps:

pt = Table[
Plot[ff[x, t], {x, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}}], {t, 0, 4, 2}];

Show[pt]


The steady state solution g[x] is given by a solution of the ODE

eqs = x g[x] + g[x] g'[x] == const

(*
Out[23]= x g[x] + g[x] Derivative[1][g][x] == const
*)


The constant can be calculated from the boundary conditions or from the condition of symmetry g'[0]==0 to give const = 0. Hence we find g[x] = 0 or, applying the boundary condition,

g[x_] := (1 - x^2)/2


The transisiton from the initial state to the steady state requires some discussion, as it is not as smooth as it might seem from the plots.

Indeed, a 3D plot with "continuous" time steps shows some instability at t ~= 0.15

Plot3D[ff[x, t], {t, 0, 1}, {x, -1, 1}]


Notice that this is the region in t where the solution makes the transisiton from positive second spatial derivative to the negative one of the steady state solution. This will be investiged later in more detail.

Conservation law

An interesting quantity is the integral w over the solution over the complete x-interval. This could be interpreted as the total thermal energy in the system

w [t_] := Integrate[f[x, t], {x, -1, 1}]


From the original PDE it might seem that this quantity must not change in time. Indeed, integrating the PDE over x from -1 to +1 and changing the order of the temporal derivative and the spatial integral gives

D[w, t] = q[x, t]/x -> +1 - q[x, t]/x -> -1


where

q[x_, t_] = x f[x, t] + f[x, t] D[f[x, t], x];


But because ot the boundary conditions, q is zero.

On the other hand w assumes different values for t = 0 and t->[Infinity], so it must have changed in the meantime. This strange behaviour will also be examined in more detail, of course.

• Your final plot is much better behaved with MaxStepSize -> 0.06. However, your "conserved quantity" is not, indicating that NDSolve is not producing a sufficiently accurate solution. Unfortunately, decreasing MaxStepSize -> 0.05 causes NDSolve to fail at t == 0.1478. Nice analysis, anyway (+1). Apr 19, 2016 at 12:43
• It seems to me that the PDE described in Edit #2 differs from that in the question in two important ways. First, because the PDE in the question is nonlinear, it can evolve to a singularity in finite t, but the linear PDE in Edit #2 cannot. Second, the conservation law derived in your original answer no longer holds for the Edit #2 PDE. Apr 23, 2016 at 19:08

As I noted in a comment above, it is quite possible for the solution of a nonlinear PDE to become singular at finite t, and that appears to be occurring here. And, as noted by Dr. Wolfgang Hintze, the right side of the PDE, when integrated over {x, -6, 6} is zero. So, the integral of f must be a constant, and indeed it is. For the parameters given in the question,

Integrate[sol[x, 0], {x, -6, 6}]
(* 2.50663 *)
ListLinePlot[Table[Integrate[sol[x, t], {x, -6, 6}], {t, 0, time, .1}],
DataRange -> {0, time}, PlotRange -> {Automatic, {2.5065, 2.5067}}]


The tiny departure from a constant value is, presumably, due to numerical inaccuracy as t approaches time, the approximate location of the singularity.

On the other hand, the steady-state solution mentioned in the question and some of the comments is 18 - x^2/2, which when integrated over {x, -6, 6} is 144. So, this steady state solution is inaccessible from the initial conditions given in the question. I should add that merely selecting initial conditions that integrate to 144 is not necessarily sufficient to reach the steady state either.

• It was just because of the "energy leak paradox" that I was absent for some days. It was new to me, as it is absent in the linear diffusion problem. But I have solved it now. Please see my EDIT. Apr 22, 2016 at 9:14