5
$\begingroup$

I see lots of cases where DSolve gives solution that has $e^{c_1}$ in the solution, where $c_1$ is the integration constant. When solving this by hand, we normally replace $e^{c_1}$ by new constant. I always wondered why DSolve does not do this automatically?

Here is an example

ode = 2*x + 3*y[x] + (3*x + 2*y[x])*y'[x] == 0;
sol = DSolve[ode, y[x], x]

Mathematica graphics

This could be made much simpler using

 sol = sol //. 4 E^C[1] -> C[1]

Mathematica graphics

Both solutions are of course correct. But my question is why Mathematica does not do this simplification automatically. It seems like something that the DSolve code could do easily using a pattern at the end, or at the moment such terms show.

Update

This is comment to the good Answer by Akku14 below. I just wanted to say that at the end, it should not make difference if one keeps the constant as $e^{c_1}$ or replace it by a new constant $c_2$. To show this, using the ode given with some i.c., this below shows that the final IVP should come out the same (too small to write in comment to the above answer, so adding it here).

Solve $$ 2x+3y+\left( 3x+2y\right) y^{\prime}=0 $$ With initial conditions $$ y\left( 0\right) =0 $$ Mathematica gives the general solution as (Using one of the two solutions it gives) \begin{equation} y=\frac{1}{2}\left( -3x-\sqrt{4e^{c_{1}}+5x^{2}}\right) \tag{1} \end{equation} Let say we keep $e^{c_{1}}$ as is. Applying initial conditions gives the equation to solve for $c_{1}$ $$ 0=-\sqrt{4e^{c_{1}}} $$ Actually Solve here did not solve for $c_1$ when asked. Even for reals. But we know that the solution must be $c_{1}=-\infty$. This is a valid constant of integration.

Substituting this into (1) gives the IVP solution as \begin{align} y & =\frac{1}{2}\left( -3x-\sqrt{4e^{-\infty}+5x^{2}}\right) \nonumber\\ & =\frac{1}{2}\left( -3x-\sqrt{5x^{2}}\right) \tag{2} \end{align} Now lets say we replaced $e^{c_{1}}$ by new constant, say $c_{2}$ before solving for the constant, which is what we would do when solving this by hand. So (1) now becomes \begin{equation} y=\frac{1}{2}\left( -3x-\sqrt{4c_{2}+5x^{2}}\right) \tag{3} \end{equation} Now applying initial conditions to (3) gives $$ 0=-\sqrt{4c_{2}} $$ This has the solution $c_{2}=0$. Substituting this in (3) gives \begin{equation} y=\frac{1}{2}\left( -3x-\sqrt{5x^{2}}\right) \tag{4} \end{equation} We see that (4) and (2) are the same IVP solution. So the only difference is the value of the constant of integration found, but this should not affect the final IVP solution.

Actually, when keeping the constant as Mathematica shows it, which is $e^{c_1}$ and solving the ODE with with these initial conditions, it gives warning:

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

But it also did give the above IVP solution at the end. But with the above warning.

But if Mathematica had changed $e^{c_1}$ to $c_1$ and only after that solved for the constant of integration, then such a warning would not even show up.

May be this is another good reason for this change in how it writes the solution?

$\endgroup$
4
  • 1
    $\begingroup$ I suppose one might extend f[C[1]] -> C[1] to any function f. I will add that 4 E^C[1] -> C[1] is not as safe as C[1] -> Log[C[1]/4] (e.g., should C[1] happen to appear in two places in different expressions). $\endgroup$
    – Michael E2
    May 30, 2022 at 15:03
  • 3
    $\begingroup$ Here's another approach: Simplify[sol, TransformationFunctions -> {Automatic, # /. c_C :> Log[c] &, # /. Cases[#, a_?NumericQ c_C :> (c :> c/a), Infinity] &, # /. Cases[#, a_?NumericQ + c_C :> (c :> c - a), Infinity] &}]. It works on DSolve[x y''[x] == y'[x] y[x], y[x], x], too. $\endgroup$
    – Michael E2
    May 30, 2022 at 15:21
  • 1
    $\begingroup$ @MichaelE2 This is nice transformation. May be you could post it as answer. DSolve should use it :) $\endgroup$
    – Nasser
    May 30, 2022 at 15:35
  • $\begingroup$ Be careful of function ranges. $\endgroup$
    – TLW
    May 31, 2022 at 1:04

2 Answers 2

5
$\begingroup$

I wasn't able quickly to find an example in which the generalization of the sort of simplification the OP wants is problematic. I recall an instance in which the constant parameters appeared in one place as something like C[1]^2 C[2] as well as in other places. Since C[1]^2 C[2] is a constant, might it be replaced by C[1] or C[2]? Not necessarily. If it's the only place C[2] appears, then you probably cannot replace it by C[1] and get the general solution. Here is a contrived case in which replacement is slightly problematic:

DSolve[x y'[x] == -x ProductLog[E^(y[x]/x)/x] + y[x], y[x], x]
(*  {{y[x] -> E^C[1] + x C[1]}}  *)

I suppose there's nothing wrong with replacing C[1] by Log[C[1]] here, but I'd prefer Exp[C[1]] to Log[C[1]]. In any case, it starts to look to me that a general, programmatic approach to simplifying constants will be a headache.

That said, there are some transformations that are probably worth trying:

# /. c_C :> Log[c] &, (* simplify E^C[1] *)
# /.                  (* simplify a * C[1] *)
  Cases[#, a_?NumericQ c_C /; a != 0 :> (c :> c/a), Infinity] &,
# /.                  (* simplify a + C[1] *)
  Cases[#, a_?NumericQ + c_C :> (c :> c - a), Infinity] &,

Below is a function that also simplifies solutions returned in the form of Function[..]. Update: The original version used Simplify, but Simplify can operate on parts of an expression. That means that one transformation might be done on one part but not the rest. This could yield an invalid solution. The code below applies each transformation to the whole expression. Each transformation plus the identity transformation is applied at each step, and the result with the fewest leaves is selected. This is repeated until no further change occurs. Update 2: I should have mentioned that rest___ represent assumptions and options that might be passed to Simplify. Also, I believe the # in the fifth (the no-substitution alternative) does not need to be simplified, since the expression will be simplified in the other three cases even if none of the substitutions are applied. If there's an edge that proves me wrong, I haven't yet thought of it. Note that the multiple Simplify[] calls will not usually cost much time since Simplify[] caches its results.

constSimplify // ClearAll;
constSimplify[dsol_, rest___] := Activate[
   FixedPoint[
    #[[First@OrderingBy[#, LeafCount, 1]]] &[
      {#,
       Simplify[# /. c_C :> Log[c], rest],
       Simplify[# /. 
         Cases[#, a_?NumericQ c_C /; a != 0 :> (c :> c/a), Infinity], 
        rest],
       Simplify[# /. 
         Cases[#, a_?NumericQ + c_C :> (c :> c - a), Infinity], 
        rest]}
      ] &,
    Inactivate[dsol, Function],
    100], Function];

Examples:

ode = 2*x + 3*y[x] + (3*x + 2*y[x])*y'[x] == 0;
sol = DSolve[ode, y[x], x]
constSimplify[sol]
sol = DSolve[x y''[x] == y'[x] y[x], y, x]
constSimplify[sol]

Mathematica graphics

Notes:

  1. It is assumed no localization of the variables in Function is needed because the solution is assumed to come from DSolve, which does not localize variables.

  2. It is not documented but it appears that OrderingBy[expr, func, 1] returns the first position of the least element (as measured by func). On this assumption, the sequence of expressions produced by the FixedPoint iterations will strictly decrease in LeafCount until no transformation produces an expression with fewer leaves than the current one (unless 100 iterations are reached).

  3. The OP's example has two independent solutions. They could be simplified independently, but constSimplify[sol] does not do that. However, constSimplify /@ sol would.


Update 2: Alternative

In response to @Nasser's comment and request for a constSimplify[] that does not simplify if no constant simplification is applied, below is a modification, constSimplify2[]. It's a little trickier, since for the given example sol = DSolve[x^3*y'''[x] + 6*x^2*y''[x] + 7*x*y'[x] + y[x] == 0, y[x], x], Simplify[sol /. c_C :> Log[c]] has a smaller leaf count than the unsimplified sol; but the logarithms do not get simplified away, and it looks horrible. Based the leaf count, the replacement seems a success and would be returned, to the chagrin of the user. To address this, I put a restriction on the application of Log to C[], namely that _C has to be in an exponent. It's not foolproof, since E^expr means any C[] in expr is "logged." Instead of the strategy mentioned in a comment below, I use Replace to return the original expression if the rule is not applied. This looks a bit more complicated, but it achieves the requested behavior.

constSimplify2 // ClearAll;
constSimplify2[dsol_, rest___] := 
  Activate[FixedPoint[#[[First@OrderingBy[#, LeafCount, 1]]] &[{
       #,
       Replace[#,
        s_ /; ! FreeQ[s, Power[_, p_ /; ! FreeQ[p, _C]]] :>
         Simplify[# /. 
           Cases[#, 
            Power[_, 
              p_ /; ! FreeQ[p, _C]] :> (c : 
                Alternatives @@ Cases[p, _C, {0, Infinity}] :> 
               Log[c]), Infinity], rest]
        ],
       Replace[#,
        s_ /; ! FreeQ[s, a_?NumericQ c_C] :>
         Simplify[# /. 
           Cases[#, a_?NumericQ c_C /; a != 0 :> (c :> c/a), 
            Infinity], rest]
        ],
       Replace[#,
        s_ /; ! FreeQ[s, a_?NumericQ + c_C] :>
         Simplify[# /. 
           Cases[#, a_?NumericQ + c_C :> (c :> c - a), Infinity], rest]
        ]
       }] &,
    Inactivate[dsol, Function], 100], Function];
$\endgroup$
11
  • $\begingroup$ I have been using your constSimplify but still do not fully understand how it works. I was wondering how hard it is to make keep the original solution as is, if there is no transformation found or needed? Here is an example of what I mean. sol=DSolve[x^3*y'''[x]+6*x^2*y''[x]+7*x*y'[x]+y[x]==0,y[x],x] whiich gives 3 terms in the solution. There is no transformation needed. But when calling constSimplify[sol] the solution is return as the original ofcourse, but the form has changed. I'd like to keep the original form, that is all. Here is screen shot also... $\endgroup$
    – Nasser
    Jul 22, 2022 at 8:25
  • $\begingroup$ !Mathematica graphics I wonder if your function can at the end detected if it did any transformations, and if not, then simply return the original solution as is? May be have a flag that is set to true if it finds it needs to do a transformation on the constant of the integration as it run? It is not a big problem, it just will be nice to keep the original solution as is, if no transformation is done. (I tried to do it myself, but I still do not understand fully your code, even after reading your notes. little too advanced for me :) Thanks. $\endgroup$
    – Nasser
    Jul 22, 2022 at 8:27
  • $\begingroup$ @Nasser It might be hard. It tries various transformations, including no transformation with no simplification. Then it picks the one with the fewest leaves (like Simplify does), which may or may not be a transformed one. I can keep track if a transformation was actually tried, but it's much harder to tell if that transformation was chosen. One could add a check at the end, e.g. Simplify[input - ouput] == 0 and return input if zero, maybe TimeConstrained. (AFAIR, using Simplify + TransformationsFunctions did not work reliably, so I opted for FixedPoint + OrderingBy.) $\endgroup$
    – Michael E2
    Jul 22, 2022 at 13:57
  • 1
    $\begingroup$ @Nasser I had an idea. I could keep a flag associated with each result indicating whether a transformation was performed. The update an overall flag with each iteration. At the end, we'll know if any substitutions succeeded. I have to go soon, so it'll take me some time to think and work on it. $\endgroup$
    – Michael E2
    Jul 22, 2022 at 14:07
  • 1
    $\begingroup$ @Nasser I came up with an alternative. I don't have a lot of test cases at hand, so let me know if it fails for you. $\endgroup$
    – Michael E2
    Jul 22, 2022 at 23:01
1
$\begingroup$

With the solutions, DSolve gives, you can access all real solutions for real x with real c1.

Substituting E^c1 -> c2 and restricting c2 also to real values, you will get some range of x, where solution is no more real.

ode = 2*x + 3*y[x] + (3*x + 2*y[x])*y'[x] == 0;

ysol1[c1_] = y /. DSolve[ode, y, x] /. C[1] -> c1

ysol2[c2_] = y /. DSolve[ode, y, x] /. E^C[1] -> c2

red1 = Reduce[
  Element[x, Reals] && Element[c1, Reals] && 
   Element[Through[ysol1[c1][x]], Reals], x]

(*   c1 \[Element] Reals && x \[Element] Reals   *)

red2 = Reduce[
  Element[x, Reals] && Element[c2, Reals] && 
   Element[Through[ysol2[c2][x]], Reals], x, Backsubstitution -> True]

(*   (x \[Element] Reals && c2 <= 0 && 
   x <= -((2 Sqrt[-c2])/Sqrt[5])) || (x \[Element] Reals && c2 <= 0 &&
    x >= (2 Sqrt[-c2])/Sqrt[5]) || (x \[Element] Reals && c2 > 0)   *)

Manipulate[
 Plot[Evaluate@Through[ysol1[c1][x]], {x, -100, 100}, 
  PlotRange -> {{-100, 100}, {-100, 100}}], {{c1, -3}, -10, 10, 
  Appearance -> "Labeled"}]

Manipulate[
 Plot[Evaluate@Through[ysol2[c2][x]], {x, -50, 50}, 
  PlotRange -> {{-50, 50}, {-50, 50}}], {{c2, -25}, -50, 50, 
  Appearance -> "Labeled"}]

enter image description here

Edit

What i really want to say: If you use solution ysol[c2][x], people may like to see, how solution behaves with different Real c2, (like my Manipulate[..ysol2[...]) and would be astonished to see, function is not real for all x. The reason is, while using negative real c2, it is as using imaginary c1.

Reduce[{E^c1 == c2, c2 == -25}, c1] 

(*   C[1] \[Element] Integers && c2 == -25 &&   c1 == I \[Pi] + 2 I \[Pi] C[1] + 2 Log[5]   *)
$\endgroup$
2
  • $\begingroup$ Thanks for the answer. You bring good point, but I think at the end of the day, the final IVP solution should be the same. I hope you agree on this. Too small to give an example here, so I added it at the end of my question. I used the ode you used and just added some initial conditions to show this. $\endgroup$
    – Nasser
    Jul 23, 2022 at 8:19
  • $\begingroup$ @Nasser , yes, of course the IVP yields the same. But you only have c1 and c2, if you don't solve the IVP. What i really wannted to say: If you use solution ysol[c2][x], people may like to see, how solution behaves with different Real c2, (like my Manipulate[.....]) and would be astonished to see, function is not real for all x. The reason is, while using negative real c2, it is as using imaginary c1. Reduce[{E^c1 == c2, c2 == -25}, c1] yields C[1] \[Element] Integers && c2 == -25 && c1 == I \[Pi] + 2 I \[Pi] C[1] + 2 Log[5] $\endgroup$
    – Akku14
    Jul 23, 2022 at 13:13

Your Answer

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

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