Bizarre solutions of a first order PDE

Question

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?

Answer 1

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]

Answer 2

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 :()