Why this PDE causes internal 1/0 division?

Question

Bug introduced in 11.3 and persisting through 12.0.0
In 11.2 DSolve return unevaluated.
Assigned WRI case number 4210941

---

Using 11.3 on Windows 10, this input

ClearAll[w, x, y, a, b, lambda, mu];

pde = D[w[x, y], x] + (a*Exp[lambda*x]*y^2 + mu*y - 
        a*b^2*Exp[(lambda + 2 mu)*x])*D[w[x, y], y] == 0;

sol = DSolve[pde, w[x, y], {x, y}]

gives

When changing 2 mu to just mu the kernel error goes away

ClearAll[w, x, y, a, b, lambda, mu];
pde = D[w[x, y], x] + (a*Exp[lambda*x]*y^2 + mu*y - 
       a*b^2*Exp[(lambda + mu)*x])*D[w[x, y], y] == 0;

sol = DSolve[pde, w[x, y], {x, y}]

why does DSolve generate this kernel error on first input above but not the second? Is this a bug?

Answer 1

DSolve undoubtedly is attempting to use the usual Method of Characteristics to solve pde. This would begin as follows:

y'[x] == (Coefficient[pde // First, D[w[x, y], y]] /. y -> y[x])
(* y'[x] == -a b^2 E^((lambda + 2 mu) x) + mu y[x] + a E^(lambda x) y[x]^2 *)
DSolve[%, y[x], x]

Unfortunately, DSolve again fails with many of the same error messages as for the PDE, returning

(* {{y[x] -> Indeterminate}} *)

Certainly, the inability of DSolve to produce a solution for any particular differential equation should not be considered a bug. However, DSolve should be trapping internal errors and should return unevaluated instead of returning y[x] -> Indeterminate. I would agree with Nasser that this is a bug, although a minor one.

Incidentally, DSolve can obtain a solution for

sol = DSolve[pde /. {mu -> 1, lambda -> 1}, w[x, y], {x, y}] // Flatten
(* {w[x, y] -> C[1][-(1/2) I (a b E^(2 x) + 2 ArcTanh[(E^-x y)/b])]} *)

and perhaps for other particular values of mu and lambda.

Addendum: General Solution Obtained Empirically

In fact, DSolve solves pde for all rational values I have tried for mu and lambda. From comparing a few such solutions, a general solution can be guessed.

(* {w[x, y] -> C[1][-I*(a*b*E^((mu + lambda) x)/(mu + lambda) + 
    ArcTanh[y E^(-mu x)/b])]]} *)

which can be verified by substituting this solution into pde. Irrational numbers seem to work too, provided that they do not have too many digits after the decimal point.

Addendum: General Solution Derived by Method of Characteristics

As shown in the fist code block of this answer, DSolve fails to solve the ODE for characteristics of pde. However, with the transformation y -> Sqrt[-b^2]*E^(mu x)*z[x], progress can be made.

ode = y'[x] == -a b^2 E^((lambda + 2 mu) x) + mu y[x] + a E^(lambda x) y[x]^2
Simplify[% /. y -> Function[{x}, Sqrt[-b^2]*E^(mu x)*z[x]], b > 0]
(* a b E^((lambda + 2 mu) x) + a b E^((lambda + 2 mu) x) z[x]^2 + 
   I E^(mu x) z'[x] == 0 *)

(Because b enters the equation only quadratically, it can be assumed positive without loss of generality. Next, apply DSolve and transform from z back to y.

Flatten@DSolve[%, z[x], x] /. Rule -> Equal /. z[x] -> y[x]/(I b*E^(mu x))
(* {-((I E^(-mu x) y[x])/b) == I Tanh[(a b E^((lambda + mu) x))/(lambda + mu)
    - I C[1]]} *)

The Method of Characteristics requires C[1] in terms of the variables.

Simplify[Solve[%, C[1]], C[2] \[Element] Integers] // Flatten
(* {C[1] -> -((I a b E^((lambda + mu) x))/(lambda + mu)) - 
    I ArcTanh[(E^(-mu x) y[x])/b] - Pi C[2]} *)

Because C[1] is an arbitrary constant, Pi C[2] can be absorbed into it (or, equivalently set to zero).

{C[1] -> -I ((a b E^((lambda + mu) x))/(lambda + mu) + ArcTanh[(E^(-mu x) y[x])/b])}

According to the Method of Characteristics, w is an arbitrary function of the right side of the expression immediately above, which is equivalent to the expression given in the first addendum. (Note that C[1] has different meanings in the two addenda, because different equations are being solved - pde and the corresponding ode.)