# Numerical solution to a nonlinear Ordinary Differential Equation

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:= eq[-0.9]

During evaluation of In:= 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= {}


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

• The error message suggests that you give "initial conditions for both values and derivatives of the functions". As far as I can see, you have an initial condition on the value of the function. Would you be in a position to give one on the derivative as well, and retry? – MarcoB Aug 1 '17 at 16:16
• I did the following: I gave an initial condition on $F'$ and obtained that for $F$ from the equation. Even so, I get error messages NDSolve::ivcon: The given initial conditions were not consistent with the differential-algebraic equations. NDSolve will attempt to correct the values. – Oliver Fabio Piattella Aug 1 '17 at 16:57

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", "F'"},
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 for one of them, and the corresponding value of F 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. == -F *)


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 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')^2 (1 - 1/(1 + F')^0.9)^2 == -F *)
F'[x] /. FindRoot[0 == sx + x /. gamma -> -0.9, {F'[x], -.8, -.9}];
f01 = -First[%%] /. F' -> %
(* -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.

• @OliverFabioPiattella The function in FindRoot is equivalent to ((x + sx) // Simplify) /. F'[x] -> fp. I can modify the answer to make this more transparent, if you like. With respect to your question about the first nonlinear solution, F[x] is equal to x/2 at small x, as can be seen from the third figure. So, F[x] is equal to x^2/4` there. – bbgodfrey Aug 2 '17 at 17:48
• Thanks a lot, it's all perfectly clear. I just needed to do some calculations (on the paper) to verify all. – Oliver Fabio Piattella Aug 2 '17 at 23:10