Issue reported to Wolfram, Inc as a possible bug in Version 12.1.1; CASE:4630268.

With Mathematica "12.1.1 for Microsoft Windows (64-bit) (June 19, 2020)", DSolve produces unexpected error messages:

 DSolve[{p1'[x] == p1[x]^2 + 2 p1[x] p2[x], 
         p2'[x] == 2 p1[x] p2[x] + p2[x]^2}, {p1, p2}, x]

Union::normal: Nonatomic expression expected at position 2 in { ... }⋃$Failed.

Flatten::normal: Nonatomic expression expected at position 1 in Flatten[$Failed].

and returns unevaluated after a few minutes. Evidently, DSolve has passed the bad argument {...}⋃$Failed to Union. I am asking

Update: Issue resolved in Mathematica "12.3.1 for Microsoft Windows (64-bit) (June 19, 2021)"

This bug has been eliminated in the sense that the internal errors no longer are leaked. Instead, DSolve simply returns unevaluated after a few minutes.

  • 2
    $\begingroup$ With v12.1.1 on a Mac I see two internal errors Union::normal and Flatten::normal $\endgroup$
    – Bob Hanlon
    Sep 9, 2020 at 13:14
  • 1
    $\begingroup$ @BobHanlon Upon further consideration, I added the second error. Thanks for the recommendation. $\endgroup$
    – bbgodfrey
    Sep 11, 2020 at 3:52

1 Answer 1


Does this count as a workaround?:

(* dividing 2nd ODE by 1st yields a homogeneous ODE *)
p2sol = DSolve[
     {p2'[p1] == (2 p1 p2[p1] + p2[p1]^2) /
                  (p1^2 + 2 p1 p2[p1])}, p2, p1] /.
    C[1] -> Log[C[1]] /. 
   p_Power :> RuleCondition[p, True];

(* p2sol turns the p1'[x] ODE in the system into a separable equation *)
PrintTemporary@Dynamic[foo = Clock[Infinity]];
   (Print[foo]; #) &@
    Flatten@ DSolve[#, p2, x],
   Print[Style[foo, Red]]; $Failed] & /@
 (Last[system] /. p1[x] -> p1[p2[x]] /. p2sol)
{{p1 -> Function[{x}, 
         Integrate][(-9 C[1] K[1]^2 + 
          Sqrt[3] Sqrt[C[1]^2 K[1]^3 (4 C[1] + 27 K[1])])^(1/3)/(
        K[1] (-2 2^(1/3) 3^(2/3) C[1] K[1] + 
           9 K[1] (-9 C[1] K[1]^2 + 
              Sqrt[3] Sqrt[C[1]^2 K[1]^3 (4 C[1] + 27 K[1])])^(1/3) + 
           2^(2/3) 3^(
            1/3) (-9 C[1] K[1]^2 + 
              Sqrt[3] Sqrt[C[1]^2 K[1]^3 (4 C[1] + 27 K[1])])^(
            2/3))), {K[1], 1, #1}] &][
     x/3 + C[2]]]}, $Failed, $Failed}

One can combine with with p2sol to obtain p2'[x]. The solutions can be stated as implicit equations, but Mathematica tries really hard to solve them.

Note this system and the one it came from admit two one-parameter families of symmetries, scalings {p1, p2, 1/x} -> C[1] {p1, p2, 1/x} and translations x -> x + C[2]. Thus they can theoretically be expressed as successive quadratures as above, provided one can solve the intermediate equations such as the one produced by this generalization:

  {p2'[p1] == (2 a p1 p2[p1] + b p2[p1]^2)/(c p1^2 + 2 d p1 p2[p1])}, 
  p2, p1]

Solve[(-c (b - 2 d) Log[p2[p1]/p1] +
   (b c - 4 a d) Log[-2 a + c - (b p2[p1])/p1 + (2 d p2[p1])/ p1]) / 
  ((2 a - c) (b - 2 d)) == C[1] - Log[p1], p2[p1]]
  • $\begingroup$ Very nice answer (+1), A few questions: x -> x + C[2] is essential to reducing the second order system to first order. What advantage does {p1, p2} -> C[1] {p1, p2} (actually, {p1, p2, 1/x} -> C[1] {p1, p2, 1/x}, I believe) provide? Where is RuleCondition documented, and what does /. p_Power :> RuleCondition[p, True] do in the code above? Should the undefined symbol p1sol actually be p2sol? $\endgroup$
    – bbgodfrey
    Sep 11, 2020 at 15:14
  • $\begingroup$ The scaling symmetry means it's reducible to a homogeneous ODE; the translation symmetry means it's autonomous, and autonomous first-order ODEs, which you get after the first integration, are separable. RuleCondition can be found only on site; look at the oldest posts; I think @WReach explains it. I used it for in-place simplification inside Function, which is HoldAll. All it does is evaluate and replace the Power[E, Log[C[1]] that comes from the previous substitution inside Function. The power evaluates to C[1], which looks much nicer. It's not important, though. $\endgroup$
    – Michael E2
    Sep 11, 2020 at 15:26
  • $\begingroup$ @bbgodfrey Yeah, p1sol should be p2sol. I solved it both ways, but it made little difference. And then I changed the fancy monitoring code and forgot to change the p1 back to p2. $\endgroup$
    – Michael E2
    Sep 11, 2020 at 15:29
  • $\begingroup$ @bbgodfrey My computer was crashing...I think my first comment was too hasty: It becomes homogeneous because of both symmetries. The system being autonomous, which is necessary to reduce the system to a single ODE, led me to ignore x or 1/x as you cleverly formulated the transformation. Scalings and homogeneity go hand-in-hand, and the single ODE obtained being homogeneous is connected to the scaling symmetry. $\endgroup$
    – Michael E2
    Sep 11, 2020 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.