I am trying to numerically solve some nonlinear partial differential equations similar to the example given below, for which I have been unable to obtain stable numerical solutions due to some peculiar features of the equations. There is a case for which a series solution exists for the steady-state solution, but I am unsure how to impose suitable boundary conditions so that the numerical solution of the full time-dependent equation is more stable and will converge to the desired steady-state solution.


An example which displays the sorts of features I have been dealing with is the following equation, which shows up in a few papers, such as this paper that I will focus on for the purposes of this post. I have modified the equation in the paper to eliminate some constants and make the time-like variable positive, and shifting their $u$ by a divergent but unimportant function of $t$ so that $\lim_{t \rightarrow +\infty} u(t,x,y) \equiv u_\ast(x,y) = 0$ is actually a valid steady-state solution, as claimed by [1]. (See the note at the end of the post for the explicit modifications I made) The resulting equation for $u(t,x,y)$ is
$$u_t(t,x,y) = (d+2) u - d~y~u_y - \frac{1 + u_{xy}}{\sqrt{(1+u_{xy})^2 - u_{xx} u_{yy}}} + 1,$$ with initial condition $u(0,x,y) \equiv u_0(x,y)$, where $u_t = \frac{\partial u}{\partial t}$, etc., and $d > 0$ is a parameter. The domains are $t \in [0,\infty)$, $y \in [0, \infty)$, and $x \in \mathbb{R}$. I should stress that although this looks a bit like a nonlinear diffusion equation, the quantity $u(t,x,y)$ is not a physical quantity like a concentration, but is a more abstract object. Notably, the IC we will look at below is unbounded as $|x|$ or $y \rightarrow \infty$.

Of particular interest are the possible fixed point solutions $u_\ast(x,y)$, and perturbations away from these solutions. Being a nonlinear equation there are many possible fixed point solutions; these are usually restricted to certain families based on symmetries of the model, or other properties of the initial condition that are preserved by the dynamics, and the fact that the fixed point solution must exist for all $y$ and $x$ on the domain. We are interested in solutions that are polynomially bounded in $y$. In particular, we expect $u_\ast(x,y) \sim y^{1+2/d}$ for large $y$. I am not sure anyone knows what the appropriate asymptotic behaviors are as $|x| \rightarrow \infty$.

The initial condition we'll focus on here is $$u_0(x,y) = g^0_{12}~x y^2/2 + g^0_{22}~x^2 y^2/4,$$ with parameters $g^0_{12},~g^0_{22} > 0.$ The authors of [1] argue that a series expansion of the solution $u(t,x,y)$ around $x = 0, y = 0$ will not contain any binomial terms $x^m y^n$ for which $m > n$. One can check that $u_\ast(x,y) = 0$ is a fixed point solution that is stable to any perturbation that obeys this restriction when $d > 2$.

When $d < 2$, the authors show the trivial solution becomes linearly unstable to a non-trivial solution, which they determine via a power series by showing that the coefficients $g_{mn} \equiv \frac{\partial^{m+n}u^\ast}{\partial x^m \partial y^n}\Big|_{x=y=0}$ obey a chain of dependencies that allows them to be evaluated iteratively: the equation for $g_{22}$ is closed, and determines $g_{12}$, which determines $g_{33}$, and so on. The pattern is $g_{n,n} \rightarrow g_{n-1,n} \rightarrow \dots \rightarrow g_{1,n} \rightarrow g_{n+1,n+1} \rightarrow \dots$, starting with $n = 2$. (See Fig. 2 of [1], and a minimal Mathematica code to compute these solutions below).

My issue: I would like to solve the PDE numerically and see the decay to $0$ for $d > 2$ and convergence to the solution that matches the series near $(x,y) = (0,0)$ for $d < 2$. A straightforward attempt at a solution does not yield good results: we use a finite range of $x$ and $y$, without explicitly imposing boundary conditions, and the solutions initially appear to be heading in the expected directions, but quickly become numerically unstable at finite $t$. Here is a minimal working example that attempts to solve the PDE:

(* PDE for u[t,x,y] with parameters d, g12, g22 *)
ucoag[d_, g12_, g22_] := 
 NDSolve[{D[u[t, x, y], t] == (d + 2)*u[t, x, y] - 
     d*y*D[u[t, x, y], y] - (1 + D[u[t, x, y], x, y])/
      Sqrt[(1 + D[u[t, x, y], x, y])^2 - 
        D[u[t, x, y], x, x]*D[u[t, x, y], y, y]] + 1, 
   u[0, x, y] == g12*x*y^2/2 + g22*x^2*y^2/4}, 
  u, {t, 0, 30}, {x, -1.0, 1.0}, {y, 0, 10}]

(* Case 1: d = 2.1 > 0. Solution expected to decay to 0. *)

ud2p1 = ucoag[2.1, 0.01, 0.01]

(* Plot result. Numerical solution shows signs of instability around t ~ 22. *)
 Plot3D[{Evaluate[u[t, x, y] /. ud2p1]}, {x, -1, 1}, {y, 0, 10}], {t, 
  0, 30}]

(* Case 2: d = 1.9 < 2. Expect convergence to non-trivial solution that matches series expansion near (x,y) = (0,0) *)

ud1p9 = ucoag[1.9, 0.01, 0.01]

(* Plot result. Numerical solution shows signs of instability around t ~ 13. *)
 Plot3D[{Evaluate[u[t, x, y] /. ud1p9]}, {x, -1, 1}, {y, 0, 10}], {t, 
  0, 30}]

(* Plot first derivative wrt x at x=0 and compare to analytical series at d = 1.9. *)

 Plot[{Evaluate[D[u[t, x, y], x] /. {x -> 0} /. ud1p9], 
   0.050000000000000044` y^2 + 0.00041666666666666767` y^3 - 
    0.000018857258812615815` y^4 + 7.675839376828031`*^-7 y^5 - 
    1.506757897829043`*^-8 y^6}, {y, 0, 0.5}], {t, 0, 12}]

Mathematica throws a warning because of the lack of boundary conditions, but otherwise will provide a solution for all $t$ specified (though the result is meaningless beyond a finite $t$ because the numerical instability will eventually cause the solution to become complex). My understanding is that Mathematica computes derivatives using only internal points when BCs are not supplied (see here and here), depending on Method. This would be similar to what has been done in solving simpler versions of similar equations numerically, using 5-point stencils.

I would assume that providing boundary conditions would help stabilize the numerical solution, but it is unclear what boundary conditions ought to be imposed. Due to the unboundedness of the solutions at one or more edges of the infinite domain, I am not sure what would be a good choice. A large constant seems too naive (and result in inconsistent boundary warnings). Perhaps implementing a polynomial bound in $y$ would help, but I am not sure how to do so.

I would appreciate any help making Mathematica's evaluation of this PDE more stable! Thank you.

Other notes:

  • [1]'s Eq. 13 reads $$u_\tau(\tau,\tilde{\chi},\chi) = -(d+2) u + d~\chi~u_\chi + \frac{\tilde{V}_d(1 + u_{\tilde{\chi}\chi})}{\sqrt{(1+u_{\tilde{\chi}\chi})^2 - u_{\tilde{\chi}\tilde{\chi}} u_{\chi\chi}}},$$ where the time-like parameter $\tau$ is negative and $\tilde{V}_d$ is a constant that depends on the parameter $d$. As written, $u(-\infty,\tilde{\chi},\chi) = 0$ is not a fixed point of this flow equation, because the authors ignore a constant offset, which does not play a role in their series solution but would diverge numerically. I explicitly removed this offset so that it does not contribute a diverging offset to the numerical solution, as well as change variables to make time positive, and to scale away the unimportant constant $\tilde{V}_d$. To be specific, to convert [1]'s Eq. 13 into the form given above and in the Mathematica code, I made the change of variables $u(\tau,\tilde{\chi},\chi) \rightarrow \tilde{V}_d u(-\tau,\tilde{\chi},\chi/\tilde{V}_d) + c(t)$, and set $t = -\tau$, $x = \tilde{\chi}$ and $y = \chi/\tilde{V}_d $. The offset can be removed by evaluating the right-hand-side at $(x,y) = (0,0)$ (for which $u(t,x=0,y=0)$ vanishes), which gives a differential equation for $c(t)$, $c_t = (d+2)c - 1$. We can then plug this into the full equation for the remaining piece $u$ to obtain the form given in this post.

  • I have tried using Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> {40, 140}}}} . I stopped evaluation after about an hour because my laptop sounded like it was going to catch fire.

  • I looked into finite element methods, but I believe this PDE cannot be put into the required standard form for the nonlinear FEM solver to work.

  • Most work solving similar equations does not actually try to solve the full PDE with $x$ dependence, instead making additional approximations or truncations to eliminate the $x$-dependence entirely. For example, similar papers define a hierarchy of equations for $v_m(t,y) \equiv \left.\frac{\partial^m}{\partial x^m}u(t,x,y) \right|_{x = 0}$, which they close via an approximation (often $v_M(t,y) = 0$ for some $M$). In this case the initial conditions become sufficient to find a solution, though there may exist artifacts due to the truncation. Hence, it would be nice to directly solve the PDE for $u(t,x,y)$, if it's possible.

  • Here is some code for iteratively calculating the coefficients $g_{mn}$. It is a bit slow for nmax larger than 6 or so, and I suspect it could be vastly improved, but that is better left as the topic of a separate question. (Unfortunately AsymptoticDSolveValue does not work with PDEs yet).

(* Determine coefficients. g[2,2] == 2(2-d), g[1,2] == 2-d determined earlier. *)

nmax = 5; (* maximum degree of x and y desired in series approximation. Can take a while for larger values. *)

glist = DeleteCases[
   Table[D[0 == ((d + 2)*u[x, y] - 
             d*y*D[u[x, y], y] - (1 + D[u[x, y], x, y])/
              Sqrt[(1 + D[u[x, y], x, y])^2 - 
                D[u[x, y], x, x]*D[u[x, y], y, y]] + 1) /. {u -> 
            Function[{x, y}, 
             Sum[g[k, l]*x^k*y^l/k!/l!, {l, 2, nmax}, {k, 1, 
               l}]]}, {x, m}, {y, n}] /. {x -> 0, y -> 0}, {n, 1, 
        nmax}, {m, 1, n}] /. {g[2, 2] -> -2 (-2 + d), 
       g[1, 2] -> -(-2 + d)} // Simplify // Flatten, True];

(* Solve for the remaining coefficients *)
gsoln = Solve[glist, 
      Table[g[m, n], {n, 2, nmax}, {m, 1, n}], {g[1, 2], g[2, 2]}] // 
     Flatten] // Flatten;

(* Plot the solution to the chosen order for, e.g., d = 1.9 *)
Plot3D[{Sum[(g[k, l] /. {g[1, 2] -> 2 - d , g[2, 2] -> 2 (2 - d)} /. 
       gsoln /. {d -> 1.9})*x^k*y^l/k!/l!, {l, 2, nmax}, {k, 1, 
    l}]}, {x, -1, 0}, {y, 0, 3}]
  • 1
    $\begingroup$ There is open access to the paper linked on arxiv.org/pdf/1302.0481v1.pdf $\endgroup$ Commented Apr 26, 2022 at 3:40
  • $\begingroup$ In the paper they proposed eq (13) that is different from what you shown above. $\endgroup$ Commented Apr 26, 2022 at 6:59
  • $\begingroup$ Yes, as I said, "I have modified the equation in the paper to eliminate some constants and make the time-like variable positive, and shifting their u by a divergent but unimportant function of t so that $\lim_{t\rightarrow \infty} u(t,x,y) \equiv u_\ast(x,y)=0$ is actually a valid steady-state solution." The 1st set of changes are only cosmetic, while the last fixes their equation to make 0 a steady-state solution, which they claim to be the case but isn't in Eq. 13. The missing 1 in their Eq doesn't affect their power series analysis, but would lead to a diverging constant offset numerically $\endgroup$
    – bbrink
    Commented Apr 26, 2022 at 11:28
  • 1
    $\begingroup$ It is not clear why in NDSolve you use equation in the form D[u[t, x, y], t] == (d + 2)*u[t, x, y]-... while in the Latex form $u_t=(d-2)u-...$? Where is typo? $\endgroup$ Commented Apr 26, 2022 at 12:46
  • 1
    $\begingroup$ The time-like variable in the paper is negative, defined to run from $0$ to $-\infty$. I defined it to be positive in this post and the Mathematica code. For the modified equation with $u_\ast = 0$ as a steady state solution, a linear perturbation $u = 0 + \delta u(t,x,y)$ yields $\delta u(t,x,y) = e^{t(d+2)} \delta u(0,x,y e^{-dt})$ for $t \geq 0$. The initial condition is quadratic in $y$, & one can show self-consistently it will not generate any linear $y$ terms, so an analytic perturbation will grow as $e^{t(d+2)} y^2 e^{-2dt} \sim e^{t(2-d)}$ at lowest order, which decays for $d > 2$. $\endgroup$
    – bbrink
    Commented Apr 26, 2022 at 16:15

2 Answers 2


We can use properties of asymptotic solution $u_*$ and that solution of equation $u_t=(d+2)u-d y u_y+V_i(t,x,y)$ defined well with polynomial initial condition and boundary condition u[t,x,0]==0. Then we have a code with iterative algorithm implementation (as Ulrich Neumann proposed) in the form


d = 3; u0 = g12*x*y^2/2 + g22*x^2*y^2/4; 
V[0] = (-(1 + D[u0, x, y])/
    Sqrt[(1 + D[u0, x, y])^2 - D[u0, x, x]*D[u0, y, y]] + 1);

g12 = g22 = .01; tmax = 5; nn = 3;

Do[sol[i] = 
  NDSolve[{D[u[t, x, y], t] == (d + 2)*u[t, x, y] - 
      d*y*D[u[t, x, y], y] + V[i - 1], 
    u[0, x, y] == g12*x*y^2/2 + g22*x^2*y^2/4, u[t, x, 0] == 0}, 
   u, {t, 0, tmax}, {x, -1, 1}, {y, 0, 10}]; 
 V[i] = (-(1 + D[u[t, x, y], x, y])/
      Sqrt[(1 + D[u[t, x, y], x, y])^2 - 
        D[u[t, x, y], x, x]*D[u[t, x, y], y, y]] + 1) /. 
   sol[i][[1]];, {i, 1, nn}]

Visualization of numerical solution on the last step

  Evaluate[Re[u[t, x, y] /. sol[nn][[1]]]], {x, -1, 1.0}, {y, 0, 10}, 
  PlotRange -> All, ColorFunction -> "Rainbow", PlotLabel -> t], {t, 
  0, tmax, .5}]

Figure 1

Visualization V[i] on every step at t=tmax and t=0.1

  Evaluate[Re[V[i] /. t -> tmax /. sol[i][[1]]]], {x, -1, 1.0}, {y, 0,
    10}, PlotRange -> All, ColorFunction -> "Rainbow", PlotLabel -> i,
   PlotPoints -> 50], {i, nn}]

  Evaluate[Re[V[i] /. t -> .1 /. sol[i][[1]]]], {x, -1, 1.0}, {y, 0, 
   10}, PlotRange -> All, ColorFunction -> "Rainbow", 
  PlotLabel -> i], {i, nn}]

Figure 2

  • $\begingroup$ Thank you for your attempt! It appears to work well up to tmax = 8 or so, beyond which the solution breaks down, for both $d > 2$ and $d < 2$ cases. Is there a straightforward way to improve the stability of the solution for longer times? $\endgroup$
    – bbrink
    Commented May 2, 2022 at 13:23
  • $\begingroup$ @bbrink There is a blow up problem as I understand since expression Sqrt[(1 + D[u, x, y])^2 - D[u, x, x]*D[u, y, y]] can reaches zero. $\endgroup$ Commented May 4, 2022 at 5:52

As a supplement to @AlexTrounev nice answer NestList might simplify the code of iteration:

d = 1.7; g12 = g22 = .01; tmax = 5  ; nn = 3;

sol[dd_, gg12_, gg22_] := 
NDSolveValue[{D[u[t, x, y], t] == (d + 2)*u[t, x, y] -d*y*D[u[t, x, y],y] + (-(1 + D[#[t, x, y], x, y])/Sqrt[(1 + D[#[t, x, y], x,y])^2 - D[#[t, x, y], x, x]*D[#[t, x, y], y, y]] +  1), 
u[0, x, y] == g12*x*y^2/2 + g22*x^2*y^2/4, u[t, x, 0] == 0}, 
u, {t, 0, tmax}, {x, -1, 1}, {y, 0, 10}] &

erg = NestList[sol[d, g12, g22], Function[{t, x, y}, g12*x*y^2/2 + g22*x^2*y^2/4], 5]
  • $\begingroup$ This is a nice answer (+1). What do you think about bbrink question concerning solution at t>8? $\endgroup$ Commented May 4, 2022 at 5:47
  • $\begingroup$ @AlexTrounev The blow up problem you mentioned occurs too in the "linear case" u0= NDSolveValue[{D[u[t, x, y], t] == (d + 2)*u[t, x, y] - d*y*D[u[t, x, y], y] , u[0, x, y] == g12*x*y^2/2 + g22*x^2*y^2/4, u[t, x, 0] == 0}, u, {t, 0, 15}, {x, -1, 1}, {y, 0, 10}]. I tried to find an iterative solution u~u0+du assuming du to be small and ic du[0, x, y] == 0, du[t, x, 0] == 0 but didn't succeed. $\endgroup$ Commented May 4, 2022 at 8:58

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.