# Bizarre solutions of a first order PDE

DSolve[{x y D[u[x,y],x]+(x-y)y D[u[x,y],y]+x==u[x,y],u[x,0]==x},u,{x,y}]


returns unevaluated. Removing the initial condition:

sol = DSolve[x y D[u[x, y], x] + (x - y) y D[u[x, y], y] + x == u[x, y], u, {x, y}]


gives something complicated; then

Simplify[u[x, 0] /. sol, Assumptions -> x \[Element] Reals]


returns

{\[Piecewise]   Indeterminate   x==0
-((x (-\[Infinity] Hypergeometric2F1[1,1+1/x,2+1/x,-\[Infinity]]
+ (1 + x) Hypergeometric2F1[1,1/x,1+1/x,-\[Infinity]]))/(1+x))  x>0
x + \[Infinity] C[1][-(x^2/2)]  True
}


which, if I understand it correctly, means that the solution is undefined at y=0 for real x.

I found several terms in the power series expansion "by hand", it goes like

x + xy/(1-x) + (xy)^2/((1-x)^2(1-2x)) + 2(2 - 3x)(xy)^3/((1-x)^3(1-2x)^2(1-3x))
+ 2(17 - 92x + 159x^2 - 90x^3)(xy)^4/((1-x)^4(1-2x)^3(1-3x)^2(1-4x))
+ 8(62 - 788x + 4048x^2 - 10783x^3 + 15759x^4 - 12042x^5 + 3780x^6)(xy)^5
/((1-x)^5(1-2x)^4(1-3x)^3(1-4x)^2(1-5x))


so that at least in the formal sense the solution seems to exist.

Can I somehow extract an explicit solution from all this?

• I would say not the solutions are bizarre but the PDE. The characteristics of this equation can be plotted with StreamPlot[{x*y,(x-y)*y},{x,-1,1},{y,-1,1},FrameLabel->{"x","y"}]. Some never come close to $y=0$. Where does this equation come from and why are you interested in it? Commented Jul 29, 2022 at 19:27
• @user293787 Very interesting observation! The equation is satisfied by the generating function for a recurrence from a recent mathoverflow question Commented Jul 29, 2022 at 21:20
• @user293787 and, on closer inspection, there is actually a discontinuity along y=0: oppositely directed streamlines clash there. Commented Jul 30, 2022 at 4:25

Not every PDE has an analytic solution. But you can always try a numerical solution:

reg = Disk[{0, 0}, 1];
0 reg = Rectangle[{-1, -1}, {1, 1}]
sol = u /.
NDSolve[{x y D[u[x, y], x] + (x - y) y D[u[x, y], y] + x == u[x, y],
u[x, 0] == x}, u, {x, y} \[Element] reg][[1]]
Plot3D[sol[x, y], {x, y} \[Element] reg, PlotRange -> All,
PlotPoints -> 50]


1. Your boundary condition isn't actually a boundary condition - it's merely a necessary condition. Consider the simpler case of xu'[x]==u[x]. u[0]==0 is necessary, but tells you nothing.
2. As you noticed from the comments, the Streamplot shows singular behavior for y=0, and to a lesser extent, x=0 and x=y as well.
3. Hypergeometric is literally a series expansion. Mathematica uses it when normal analytic methods fail.
4. Looking at your series expansion, I'm assuming you expanded it in powers of y. Well, your 5th term has a pole at x=1/5, and a 5th order pole at x=1. If I had to guess, you're faced with a sequence of poles that approach x=0. Maybe those terms will converge for sufficiently large x, but I doubt they're useful for real x between (0,1). It's also another hint that y=0 is highly singular.
5. I used the substitution v=u-x. This results in xy(dv/dx+1)+(x-y)y(dv/dy)=v. If it weren't for that pesky +1 term, this would be linear. Since you're interested in y=0, we can estimate local behavior in some epsilon ball of (0,0) by dumping that +1. Looking at that solution, it's clear that y=0 is highly singular.

Overall, the PDE itself is naughty, and mathematica can't fix that. Given the comments, you can give streamplot3d a try (doesn't work on wolfram alpha :()

• Thank you, this is most informative! The question though is mostly about how to understand the output of Mathematica in this case. Because it seems to produce some explicit solution but I don't know how it can be used since I cannot evaluate it anywhere. Commented Jul 30, 2022 at 13:50
• @მამუკაჯიბლაძე 1) Mathematica is known for producing illegible answers because dealing with all the poles/branch cuts/other special conditions can be very challenging. Even simpler things like sqrt and log can be surprisingly complicated. Sometimes you can simplify things by adding Assumptions->x>0 or w/e. 2) As I said, you haven't supplied enough boundary conditions, so there is no equation to evaluate. There's probably a c_1(...) floating around somewhere. Commented Jul 31, 2022 at 2:45
• I want that solution with series x+xy/(1-x)+O(y^2) that I had there towards the end of the question, but I don't know how to specify it with boundary conditions Commented Jul 31, 2022 at 5:41