First order Differential equation with boundary condition

Question

I'm wondering why it is not solving with the Mathematica??

Clear["Global`*"]
DSolve[{(y'[x])^2 == a*y[x] - (k*y[x]^2) + ((f*y[x]^4)/3), y[1] == 1}, y[x], x]

If you can find anyway to plot it for the some value for a,k and f? thank you

Answer 1

I just solved it without you boundary condition:

 Clear["Global`*"]
 eqns2 = {g'[x]^2 == a (f/3 g[x]^3 + 1 - k/a g[x]) g[x]};
 sol2 = DSolve[eqns2, g[x], x]

Answer 2

Plots can be obtained as follows. First, solve the ODE without the boundary condition.

sol = DSolve[{(y'[x])^2 == a*y[x] - (k*y[x]^2) + ((f*y[x]^4)/3)}, [x], x]

The solution is too lengthy to reproduce here but has the form,

{Solve[exp == x/Sqrt[3] + C[1], Solve[exp == -x/Sqrt[3] + C[1]}

where the lengthy expression exp is the same for both solutions. However, C[1] is not necessarily the same for both. Let us focus on the first solution. The constant of integration, C[1] is obtained from

sol[[1, 1]] /. x -> 1 /. y[1] -> 1;
c = FullSimplify[%] // First;

which also is rather lengthy. Now construct with, for instance, {a -> 1, f -> 1, k -> 1},

arg = N@SolveValues[sol[[1, 1]] /. {y[x] -> w, C[1] -> c} /. {a -> 1, f -> 1, 
     k -> 1}, x] // First // Simplify;

which again is lengthy but can be plotted without difficulty. (Note that y[x] is replaced by w for convenience.

p1 = ParametricPlot[{Chop[arg], w}, {w, 0, 4}, AxesLabel -> {x, "y[x]"}, 
    ImageSize -> Large, LabelStyle -> {15, Bold, Black}]

To obtain the second half of this curve, note that it simply is the first half, just plotted, reflected in x about the point where w = 0.

Chop[arg /. w -> 10^-12]
(* 3.48209 *)

So, the total solution is

p2 = ParametricPlot[{2 3.48209 - Chop[arg], w}, {w, 0, 4}, AxesLabel -> {x, "y[x]"}, 
    ImageSize -> Large, LabelStyle -> {15, Bold, Black}];
Show[p1, p2, PlotRange -> {{-1, 8}, Automatic}]

By the symmetry already noted of the two solutions contained in sol, the second curve passing through y[1] = 1 is identical to the plot immediately above but centered at x = -1.48209 instead of x = 3.48209.

Incidentally, it is simpler but not quite as accurate to use NDSolve to obtain one of the two branches and, once again, adding the second branch by reflecting the first about the value of x for which y[x] = 0.

Addendum

A still simpler approach for a symbolic solution is to determine x directly as a function of y, including the boundary condition.

solx = DSolve[{x'[y]^2 (a*y - k*y^2 + f*y^4/3) == 1, x[1] == 1}, 
    x[y], y] // Values // Flatten

which then can be plotted with ParametricPlot as above.

argx = solx /. {a -> 1, f -> 1, k -> 1};
ParametricPlot[Evaluate@Thread[{argx, y}], {y, 0, 4}]

Answer 3

A completely numerical solution can be obtained quite simply as follows. With

eq = y'[x]^2 == a*y[x] - k*y[x]^2 + f*y[x]^4/3

differentiate eq and eliminate common terms to obtain a second order autonomous ODE.

D[eq, x];
eqlin = DivideSides[%, 2 y'[x], Assumptions -> y'[x] != 0] // Cancel
(* (y''[x] == 1/6 (3 a - 6 k y[x] + 4 f y[x]^3) *)

The two corresponding initial conditions are

y[1] == 1

eq /. x -> 1 /. y[1] -> 1
(* y'[1]^2 == a + f/3 - k *)

If real solutions are required, as seems likely, then a + f/3 - k > 0 is necessary. For the parameters used in my earlier, symbolic solution,

Quiet@NDSolve[{eqlin, y[1] == 1, y'[1]^2 == a - k + f/3} /. {a -> 1, 
    f -> 1, k -> 1}, y[x], {x, -10, 10}] // Flatten // Values;
Join[{x}, Interval @@ Flatten[% /. x -> "Domain", 1] // First]
Quiet@Plot[%%, %, AxesLabel -> {x, "y[x]"}, ImageSize -> Large, 
    LabelStyle -> {15, Bold, Black}, PlotRange -> {0, 5}]

as expected. In general, curves are bounded from above by the smallest root greater than 1 of,

a*y - k*y^2 + f*y^4/3 == 0

if any. Otherwise, the curves are unbounded from above, as in the plot above. Similarly, the curves are bounded from below by the largest root less than 1. Note that y == 0 always is a root. Here is an example of a solution bounded from both above and below.

FindInstance[{a - k*y + f*y^3/3 == 0, y > 1, a + f/3 - k > 0}, {a, f, k, y}] 
    // N // Flatten;
param = % // Most
(* {a -> 0.887756, f -> -4.88266, k -> -2.21939} *)

Quiet@NDSolve[{eqlin, y[1] == 1, y'[1]^2 == a - k + f/3} /. param, 
    y[x], {x, -10, 10}] // Flatten // Values;
Join[{x}, Interval @@ Flatten[% /. x -> "Domain", 1] // First];
Quiet@Plot[%%, %, AxesLabel -> {x, "y[x]"}, ImageSize -> Large, 
    LabelStyle -> {15, Bold, Black}, PlotRange -> {0, 1.4}]