# Singularity while solving PDE

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?

• I get another warning, NDSolve::bcart -- why aren't you concerned with that one? Or don't you get the message? Oct 11 '15 at 12:20
• So you're solving the equation on an infinite region? Oct 11 '15 at 12:23
• This error message is quite common (unfortunately) when solving PDEs. It can have many causes, the most common of which (in my experience) is a numerical instability. Sometimes, increasing resolution helps, but often it makes things worse. Oct 11 '15 at 20:09
• @bbgodfrey In OP's case it's mainly because of the absence of the boundary condition. BTW the behavior of NDSolve when boundary conditons aren't enough isn't fully understand so far, I used to ask this but didn't get a satisfied answer. Oct 12 '15 at 2:28
• @xzczd Agreed. Too bad you did not get an answer to your question (+1). I started out with my answer below, hoping to find some simple work-around but did not. Oct 12 '15 at 2:33

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!