Message Boards Message Boards

DSolve not giving me an answer

Posted 10 months ago

I'm confused why DSolve can't solve this problem:

DSolve[{y'[x] == (7 x + 4 y[x])/(7 x + 4 y[x] + 4), y[-1] == -1}, 
 y[x], x]

Can you let me know how I can solve this using Mathematica?

POSTED BY: Henric Larsson
7 Replies

Looks like to me is a bug in DSolve. Solution is:

$y(x)=-\frac{7}{11}-\frac{7 x}{4}+\frac{4}{11} W_{-1}\left(\frac{1}{16} (-93) e^{\frac{7}{4}+\frac{121 x}{16}}\right)$

eq = {y'[x] == (7 x + 4 y[x])/(4 + 7 x + 4 y[x]), y[-1] == -1};
eq /. y -> 
   Function[{x}, -(7/11) - (7 x)/4 + 
     4/11  ProductLog[-1, -(93/16)  E^(
        7/4 + (121 x)/16)]] // FullSimplify

(*{True, True}*)

Regards M.I

POSTED BY: Mariusz Iwaniuk

The numerical solution with NDSolve misses a singular point with infinite derivative. Here is a derivation of the solution:

diffEq = y'[x] == (7  x + 4  y[x])/(7  x + 4  y[x] + 4);
solY = DSolveValue[diffEq, y, x]
eq = y == solY[x]
extractProductLog = Solve[eq, Cases[eq, _ProductLog, All][[1]]][[1, 1]]
implicitEq = 
 extractProductLog /. (ProductLog[z_] -> w_) :> z == w  Exp[w]
initialCond = implicitEq /. {x -> -1, y -> -1}
solForC1 = Solve[initialCond, C[1], Reals]
implicitEquationForInitialConditions = 
 Simplify[implicitEq /. solForC1[[1]]]
Reduce[implicitEquationForInitialConditions, y, Reals]
ContourPlot[implicitEquationForInitialConditions,
 {x, -3, 3}, {y, -3, 3}]
POSTED BY: Gianluca Gorni

The problem is that the use of inverse functions by Solve picks an invalid branch of ProductLog for the generic solution to the differential equation. One can tentatively assume that ProductLog[u] may be replaced by ProductLog[k, u], and check its validity after obtaining a solution. This yields the same result as Mariusz's.

{ode, ic} = (* I. separate ODE and IC *)
  {y'[x] == (7  x + 4  y[x])/(7  x + 4  y[x] + 4), y[-1] == -1};
DSolve[ode, y[x], x] /.(* II. general solution *)
     ProductLog -> (* III. replace principal ProductLog by general branch *)
      (ConditionalExpression[ProductLog[C[2], #], 
         C[2] \[Element] Integers] &) //
    (# /. Solve[ (* IV. solve IC for parameters *)
          ic /. DSolve`DSolveToPureFunction[#] (* Change y[x] rule to Function *)
          , {C[1], C[2]}] & /@ (* Map IC-solving over each DSolve solution *)
       # &) //
   Apply@Join //(* V. Join particular solutions *)
  FullSimplify[#, _C \[Element] Integers] & // Expand  (* VI. simplify *)

(* {{y[x] -> -7/11 - 7 x/4 + 4/11 ProductLog[-1, -93/16 E^(7/4 + 121 x/16)]}} *)

The above code is more-or-less a standard work-flow I use:

  1. If DSolve fails on IVP, separate ODE and IC, get general solution, and try for the particular solution to the IVP.

  2. Get the general solution.

  3. Replace principal-branch ProductLog by a general-branch version. For the branch number k, I used C[2]. It turns out that Solve can handle the IC with the general-branch product-log, and solve for the branch. I wonder if this could be added to how DSolve sets up the IC. It might be worth reporting to WRI.

  4. Solve the initial condition for each solution in the general solution returned by DSolve. I use Map (/@) to do this. For nonlinear ODEs, it is possible to get more than one solution. In the case at hand, there is only one, and the use of methods that avoid Map is possible. But those might fail, if the assumption that there is only one turned out to be wrong (for instance, after fixing a typo).

  5. Solvers like DSolve and Solve always yield a list of solutions. So solving the IC gives us a list of lists of solutions that need to be joined. (Again, as in point 4., there is only one in the case at hand, and there are other ways work in this case, which do not work in general.)

  6. Simplifying is optional. In the case at hand, a spurious, conditional, third integer parameter C[3] is generated by Solve (something like Log[expr] + 2 Pi I C[3], which is plugged into Exp[], making the parameter disappear). The simplifying gets rid the condition.

The steps can be inspected by removing the // and saving the result in a variable. E.g., for a // func1 // func2, write res1 = func1[a]; res2 = func2[res1]. And so forth. You can also insert // Echo to have the result printed out below the input and above the output.

POSTED BY: Michael Rogers
Posted 10 months ago

Thanks to all for your help. The solution can be found by using substitution u = 7x + 4y. Not sure if the algorithm can be improved to handle these problems.

POSTED BY: Henric Larsson
Posted 10 months ago

DSolve solves your ode

DSolve[{y'[x] == (7 x + 4 y[x])/(7 x + 4 y[x] + 4)}, y[x], x]

as an Abel equation of the second kind. The problem is with solving the initial condition for the parameter of integration. Not sure how the substitution u = 7 x + 4y helps with that.

POSTED BY: Updating Name

With substitution u = 7 x + 4y not helps.

POSTED BY: Mariusz Iwaniuk

Another way, using Henric's substitution to get an autonomous DE. The DE being first order, this reduces the solution to quadrature: $$u'=f(u) \wedge u(x_0)=u_0 \Longrightarrow \int_{u_0} {du \over f(u)} = \int_{x_0} dx \,.$$ The use of the assumption assum in a couple of steps is crucial to obtain valid and only valid solutions.

sub = y[x] == (u[x] - 7 x - 4)/4;
assum = 4 + 7  x + 4  y[x] < 0;
uode = DSolveChangeVariables[Inactive[DSolve][ode, y[x], x], u, x, sub] // Simplify
uprime = First@Solve[First[uode], u'[x]]
implsol = Assuming[Reduce[assum /. First@Solve[sub, y[x]], u[x]],
  Integrate[1/u'[x] /. uprime /. u[x] -> u
    , {u, First@SolveValues[sub /. {y[x] -> -1, x -> -1}, u[-1]], u[x]}
    , Assumptions -> $Assumptions] == Integrate[1, {x, -1, x}]
  ]
implsol = implsol /. First@Solve[sub, u[x]] // FullSimplify
Block[{Simplify = FullSimplify},
 ysols = Solve[implsol, y[x], Method -> Reduce, Assumptions -> assum]]
{ode, ic} /. DSolve`DSolveToPureFunction[ysols] //
 FullSimplify[#
   , Reduce[assum /. ysols // Normal, x, Reals]
   , TransformationFunctions -> {Automatic
     , ReplaceAll[eq_Equal :> FullSimplify[ApplySides[Exp, eq]]]}] &

Five lines of output of code, plus a warning message from Solve

The use of assum, Reduce, ApplySides, and FullSimplify above suggest the robustness required by the system to solve the problem needed some special tricks to be able to handle ProductLog in a careful and correct way. From the beginning, I thought Mathematica was meant, not to solve every problem for us, but to be a powerful computation tool to use in our problem solving. We seem to be in that realm here. However, below is a way to derive the assumption from the domain of the ODE, a somewhat natural thing to try, and, in this case, to derive the substitution:

assum = LogicalExpand@Reduce[
    FunctionDomain[y'[x] /. Solve[ode, y'[x]], {x, y[x]}] /. 
     Unequal -> Function[{x, y}, x < y || x > y],
    {x, y[x]}, Reals] //
  Replace[#, 
    intervals_Or :> Pick[intervals, List @@ intervals /. icToRules[ic, y[x]]]] &
sub = Reduce[u[x] == Subtract @@ assum, y[x]]

Also, I've known for a while that DSolve generally works by antidifferentiation, putting off the boundary conditions to the end. Antidifferentiation is easiest (or at least faster) than definite integration. So if the ode can be transformed into a form that is solvable by quadrature, as above, we can put the burden of the initial condition on Integrate, which will work hard, and sometimes long, to integrate from the condition $(x_0,y(x_0))$ to a general $(x,y(x))$, making sure the right branch cuts are followed. That is not the problem here. The problem is not in integrating the ODE; rather, it's in solving the implicit solution. In an intro diff. eq. class, we usually accept the implicit solution, when solving it requires special functions. In our case at hand, we need not use Integrate. We can use DSolve. We can get the same implsol as above using DSolve to solve for x[u] instead of for u[x]:

xode = First@uode /. {u'[x] -> 1/x'[u], u[x] -> u, x -> x[u]}
implsol = First@First@DSolve[
      {xode, x[-7] == -1}, x[u], u,
      Assumptions -> assum /. First@Solve[sub, y[x]] /. u[x] -> u
      ] /. {x[u] -> x, u -> u[x], Rule -> Equal};
implsol = implsol /. First@Solve[sub, u[x]] // FullSimplify
POSTED BY: Michael Rogers
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract