盐城射阳长荡镇把“雷锋火种”送到敬老院

Question

I am trying to plot the solutions to a differential equation

ClearAll;
eqn = y'[x] == (2 x + y[x] + 2)/(2 x + y[x] - 4);
solution = DSolve[eqn, y[x], x]

which gives

{{y[x] -> -2 (-2 + x) - 
  2 (1 + ProductLog[-E^(-1 - (3 x)/2 + C[1])])}}

and then plot it

sol1 = Table[ySol = (y[x] /. solution) /. C[1] -> Cval, {Cval, -20, 20, 0.5}];
sol2 = Table[ySol = (y[x] /. solution) /. C[1] -> Cval - Pi*I, {Cval, -20, 20, 0.5}];
sol = Append[sol1, sol2];
Plot[sol, {x, -4, 4}, PlotRange -> {-2, 8}]

which gives

The upper integral lines stop where the derivative becomes infinite, but they should continue (bending to the right)? How can I achieve this?

Also, is there really no way to get DSolve to produce an implicit solution to this differential equation?

Answer 1

Another way: Reverse the solution to get an implicit form:

Clear[c];
impl = Solve[
      First@DSolve[eqn, y[x], x, GeneratedParameters -> c] /. 
       Rule -> Equal, c[1]][[1, 1]] /. Rule -> Equal // Simplify // 
  Quiet

(*  c[1] == 1 + (3 x)/2 + Log[1/2 E^(1 - x - y[x]/2) (-2 + 2 x + y[x])]  *)

ContourPlot[
 Evaluate@
  Table[impl /. Log[z_] :> Log[Abs@z] /. y[x] -> y, {c[1], -10, 10, 
    1/2}], {x, -5, 5}, {y, -2, 8}, MaxRecursion -> 3]

Answer 2

Include the excluded branch alluded to by the warning message:

sol1 = Table[(y[x] /. solution) /. C[1] -> Cval, {Cval, -20, 20, 0.5}];
sol3 = Table[(y[x] /. solution) /. C[1] -> Cval /. 
    ProductLog -> (ProductLog[-1, #] &), {Cval, -20, 20, 0.5}];
sol2 = Table[(y[x] /. solution) /. C[1] -> Cval - Pi*I, {Cval, -20, 
    20, 0.5}];
sol = Join[sol1, sol2, sol3];
Plot[sol, {x, -4, 4}, 
 PlotStyle -> Table[ColorData[97][k], {k, 2 Length@sol1}], 
 PlotRange -> {-2, 8}]

Answer 3

Equation for $y(x)$ can be converted into a set of equations for $x(t)$ and $y(t)$

eqs = {y'[t] == 1, x'[t] == (2 x[t] + y[t] - 4)/(2 x[t] + y[t] + 2)}
sols = DSolve[eqs, {y[t], x[t]}, t]

Then plot the solution

ParametricPlot[
 Table[{1/2 (-2 - t) + 2 (1 - s *ProductLog[E^(-1 + (3 t)/4 + c)]), 
    t}, {c, -30, 30, 5}, {s, -1, 1, 2}] // Flatten, {t, -30, 30}, 
 PlotRange -> {{-20, 20}, {-20, 20}}]

Answer 4

• StreamPlot

StreamPlot[{1, (2 x + y + 2)/(2 x + y - 4)}, {x, -5, 5}, {y, -5, 5}]

• ParametricNDSolveValue

Clear[sol,pts];
sol = ParametricNDSolveValue[{x'[t] == 1, 
    y'[t] == (2 x[t] + y[t] + 2)/(2 x[t] + y[t] - 4), x[0] == x0, 
    y[0] == y0}, {x[t], y[t]}, {t, 0, 100}, {x0, y0}];
pts = Join[
   RandomPoint[
    ImplicitRegion[Abs[2 x + y - 4] == 10^-2, {{x, -5, 5}, {y, -5, 5}}],
     300], RandomPoint[RegionBoundary@Rectangle[{-5, -5}, {5, 5}], 
    300]];
ParametricPlot[sol @@@ pts // Evaluate, {t, 0, 100}, PlotRange -> 5]

Answer 5

Two answers have gone the parametric route. I'd suggest separating the equations as follows, though:

eqs = {y'[t] == (2 x[t] + y[t] + 2), x'[t] == (2 x[t] + y[t] - 4)};
sols = DSolve[eqs, {y[t], x[t]}, t, GeneratedParameters -> c];
tbl = Flatten[
    Table[# /. {c[1] -> -2/3 + 2 b1 + a1, c[2] -> -2/3 - b1 - 2 a1}
     , {a1, -7, 6, 1}, {b1, -3, 3, 6}], 1] &;
ParametricPlot[
 tbl[{x[t], y[t]} /. sols[[1]]] // Evaluate
 , {t, -5, 5}, PlotRange -> {{-10, 10}, {-10, 10}},
 Epilog -> {Black, Point@tbl[2/3 + {c[1], c[2]}]}]

The epilog shows the (ahem) Point of the funny Table[] substitution in tbl.