2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – MarcoB Aug 1 '17 at 16:16
  • $\begingroup$ 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. $\endgroup$ – Oliver Fabio Piattella Aug 1 '17 at 16:57
2
$\begingroup$

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}]

enter image description here

(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]

enter image description here

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]

enter image description here

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]]

enter image description here

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]]

enter image description here

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

$\endgroup$
  • 1
    $\begingroup$ @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. $\endgroup$ – bbgodfrey Aug 2 '17 at 17:48
  • $\begingroup$ Thanks a lot, it's all perfectly clear. I just needed to do some calculations (on the paper) to verify all. $\endgroup$ – Oliver Fabio Piattella Aug 2 '17 at 23:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.