Numerical solution to a nonlinear Ordinary Differential Equation

Question

I am trying to solve the following ODE for $f(x)$:

$$x f' - f = \frac{(f')^2}{\gamma^2}[1-({f}')^\gamma]^2$$

where $\gamma < 0$ and real, which makes the ODE highly nonlinear. Since I am not looking for a linear solution ($f = x$ does solve the above ODE), I have reformulated the equation for $f(x) = x + F(x)$. Here is my input:

eq[gamma_] := 
 NDSolve[{(1 + F'[x])^2*(1 - (1 + F'[x])^gamma)^2/gamma^2 == 
    x*F'[x] - F[x], F[0.01] == 0.01}, F, {x, 0.01, 100}, 
  Method -> {"EquationSimplification" -> "Residual"}]

I had to use the suggested Method -> {"EquationSimplification" -> "Residual"} in order to solve for the derivative. When I try to get a solution:

In[75]:= eq[-0.9]

During evaluation of In[75]:= NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

Out[75]= {}

I got this error which, being honest, I don't quite fully understand. Could anyone help me?

Thanks for your attention.

Answer 1

Explanation of Error Noted in Question

The computation in the question fails, because there is no real value of F'[0.01] for F[0.01] == 0.01 and gamma == -0.9. This can be seen by solving the ODE for F[0.01] as a function of F'[0.01].

sf = f /. Flatten@Solve[((1 + F'[x])^2*(1 - (1 + F'[x])^gamma)^2/gamma^2 == 
     x*F'[x] - F[x]) /. gamma -> -0.9 /. F'[x] -> fp /. F[x] -> f /. x -> 0.01, f]
(* 0.01 fp - 1.23457 (1. + fp)^2 (1. - 1./(1. + fp)^0.9)^2 *)
FindMaximum[sf, fp]
(* {0.0000249876, {fp -> 0.00499627}} *)

So, choose F[0.01] less than that maximum, for instance,

s = ParametricNDSolveValue[{(1 + F'[x])^2*(1 - (1 + F'[x])^gamma)^2/
    gamma^2 == x*F'[x] - F[x], F[0.01] == 10^-5}, F, {x, 0.01, 100}, {gamma}, 
    Method -> {"EquationSimplification" -> "Residual"}]

Plot[s[-.9][x], {x, 0.01, 100}]

(It is convenient to use ParametricNDSolveValue for problems like this.)

General Linear Solution

Differentiating the ODE reveals that there are two categories of solutions.

der2 = Collect[Subtract @@ D[(1 + F'[x])^2*(1 - (1 + F'[x])^gamma)^2/gamma^2 == 
    x*F'[x] - F[x], x], F''[x], Simplify]
(* (-x + (2 (1 + F'[x])^(1 + gamma) (-1 + (1 + F'[x])^gamma))/gamma + 
   (2 (1 + F'[x]) (-1 + (1 + F'[x])^gamma)^2)/gamma^2) F''[x] *)

Solutions given by F''[x] == 0 are, of course, linear, F[x] == a x + b. Insert this into the ODE to obtain b as a function of a.

sb = -First[((1 + F'[x])^2*(1 - (1 + F'[x])^gamma)^2/gamma^2 == 
     x*F'[x] - F[x]) /. F[x] -> a x + b /. F'[x] -> a]
(* ((1 + a)^2 (1 - (1 + a)^gamma)^2)/gamma^2 *)

ParametricPlot[{sb /. gamma -> -0.9, a}, {a, -1, 1}, AxesLabel -> {"F[0]", "F'[0]"}, 
    ImageSize -> Large, LabelStyle -> Directive[12, Black, Bold], 
    AspectRatio -> 1/GoldenRatio]

The linear solutions are a one-dimensional family, parameterized by a or b. The solution obtained in the previous section is an example.

Nonlinear Solutions

The other factor of der2 is an algebraic expression for F'[x].

sx = First[der2]
(* -x + (2 (1 + F'[x])^(1 + gamma) (-1 + (1 + F'[x])^gamma))/gamma + 
   (2 (1 + F'[x]) (-1 + (1 + F'[x])^gamma)^2)/gamma^2 *)

ParametricPlot[{x + sx /. F'[x] -> fp /. gamma -> -0.9, fp}, {fp, -1, 1}, 
    AxesLabel -> {x, "F'[x]"}, ImageSize -> Large, 
    LabelStyle -> Directive[12, Black, Bold], AspectRatio -> 1/GoldenRatio]

For x > -1 there are two solutions. Visibly, F'[0] == 0 for one of them, and the corresponding value of F[0] is given by

((1 + F'[x])^2*(1 - (1 + F'[x])^gamma)^2/gamma^2 == 
    x*F'[x] - F[x]) /. x -> 0 /. gamma -> -0.9 /. F'[0] -> 0
(* 0. == -F[0] *)

With these values, the solution is given by

fprime[gg_?NumericQ, xx_?NumericQ] := 
    fp /. FindRoot[sx /. {F'[x] -> fp, gamma -> gg, x -> xx}, {fp, xx/2}]
Flatten@NDSolve[{F'[x] == fprime[-.9, x], F[10^-4] == 10^-8/4}, F, {x, 10^-4, 100}];
Plot[F[x] /. %, {x, 10^-4, 100}, AxesLabel -> {x, F}, ImageSize -> Large, 
    LabelStyle -> Directive[12, Black, Bold]]

F[0] is determined for the second solution by

((1 + F'[x])^2*(1 - (1 + F'[x])^gamma)^2/gamma^2 == 
    x*F'[x] - F[x]) /. x -> 0 /. gamma -> -0.9
(* 1.23457 (1 + F'[0])^2 (1 - 1/(1 + F'[0])^0.9)^2 == -F[0] *)
F'[x] /. FindRoot[0 == sx + x /. gamma -> -0.9, {F'[x], -.8, -.9}];
f01 = -First[%%] /. F'[0] -> %
(* -0.599484 *)

fprime1[gg_?NumericQ, xx_?NumericQ] := 
    fp /. FindRoot[sx /. {F'[x] -> fp, gamma -> gg, x -> xx}, {fp, -.8, -.9}]
Flatten@NDSolve[{F'[x] == fprime1[-.9, x], F[10^-4] == f01}, F, {x, 10^-4, 100}];
Plot[F[x] /. %, {x, 10^-4, 100}, AxesLabel -> {x, F}, ImageSize -> Large, 
    LabelStyle -> Directive[12, Black, Bold]]

It is nearly linear, because fprime1[-0.9, x] is nearly constant.