Singularity while solving PDE

Question

I have to solve a PDE (in the context of the functional renormalization group in physics). I have a function of two variables, $U(l,p)$. I know $U(0,p)=U_0(p)$ and then I have an equation (eq925) that gives me $\partial_l U(l,p)$. Here is the code

Definitions

eq925 = D[U[l, p], l] == d U[l, p] - (d - 2) p D[U[l, p], p] - k[d]/d (1/(1 + D[U[l, p], p] + 2 p D[U[l, p], {p, 2}]));
k[d_] := (2^(d - 1) π^(d/2) Gamma[d/2])
U0[p_] := u0/6 (p - p0)^2
u0v = 0.01;
p0v = 0.050459;
d = 3;

Attempt at solving

s2 = NDSolve[{eq925, U[0, p] == Evaluate[U0[p] /. {u0 -> u0v, p0 -> p0v}]},U[l, p], {l, 0, 1}, {p, -4, 4}]

Depending on the range of p in my NDSolve I get problems. If I only keep positive p, I get a reasonable solution, but if include (like I did in the previous line code) negative values of p, I get this error

NDSolve::eerr: Warning: scaled local spatial error estimate of
 3152.029270554484` at l = 0.0659979486053968` in the direction of independent variable 
 p is much greater than the prescribed error tolerance. Grid spacing with 25 points may 
 be too large to achieve the desired accuracy or precision. A singularity may have formed 
 or a smaller grid spacing can be specified using the MaxStepSize or MinPoints 
 method options. >>

I know I could have singularities in the denominator of eq925, but this should happen in my case for values around $p=-100$, which is far outside my range.

Any ideas why this happens?

Answer 1

This is not an answer but may provide some helpful insight. The code as written contains some unnecessary functions which slow the computation slightly and reduce clarity. Instead, one could use

d = 3; u0v = .01; p0v = 0.050459;
k = (2^(d - 1) π^(d/2) Gamma[d/2]);
U0 = u0v/6 (p - p0v)^2;
eq925 = D[U[l, p], l] == d U[l, p] - (d - 2) p D[U[l, p], p] - 
    k/d (1/(1 + D[U[l, p], p] + 2 p D[U[l, p], {p, 2}]));

s2 = First@NDSolve[{eq925, U[0, p] == U0}, U[l, p], {l, 0, 1}, {p, -1.5, 4}, 
    MaxSteps -> 100000]

-1.5 is used as the lower boundary for p, because it seems to be the approximate threshold for instability. Even so, MaxSteps -> 100000 is necessary to reach l == 1. As noted by Michael E2 in a comment, NDSolve immediately issues

NDSolve::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable p. Artificial boundary effects may be present in the solution. >>

and this is correct. The system needs two boundary conditions at constant p but has none. This may be the root of the problem. (Without the need boundary conditions, NDSolve in effect makes some up.) Later, NDSolve also issues the NDSolve::eerr: warning. To see at least the symptoms of this problem, it is helpful to plot the solution.

Plot3D[U[l, p] /. s2, {l, 0, 1}, {p, -1.5, 4}, AxesLabel -> {l, p, U}]

which looks reasonable, apart from a few blemishes. However, the denominator (1 + D[U[l, p], p] + 2 p D[U[l, p], {p, 2}]) is worrisome, so it too is plotted.

Plot3D[(1 + D[U[l, p] /. s2, p] + 2 p D[U[l, p] /. s2, {p, 2}]) /. 
    p -> q, {l, 0, 1}, {q, -1.5, 4}, AxesLabel -> {l, p, U}, PlotRange -> {-.5, 1.5}]

and it shows a developing instability. Moreover, the quantity goes through zero, first at about {0.62, -1.5}. Thereafter, the solution is meaningless. So, it is essential to provide boundary conditions, with at least one of them at the minimum value of p. I tried the obvious,

(D[U[l, p], p] /. p -> -1.5) == 0, (D[U[l, p], p] /. p -> 4)==0

but it only increased the instability.

Incidentally, replacing u0v = .01 by u0v = 1/100 also increased the instability. Bizarre!