The numerical solution of a nonlinear ODE: boundary value problem

Question

I want to solve the following nonlinear ordinary differential equation with boundary values $f(0)=1$ and $f(1)=1$:

$${\frac {{\rm d}^{2}}{{\rm d}{x}^{2}}}f \left( x \right) ={\dfrac { \left( 51-51\,x \right) \left( 1+ \left( {\dfrac {\rm d}{ {\rm d}x}}f \left( x \right) \right) ^{2} \right) ^{3/2}}{25}}+{ \frac {1+\left({\dfrac{\rm d}{{\rm d}x}}f \left( x \right)\right)^2 }{f \left( x \right) }} $$

But:

NDSolve[{y''[x] == 
   204/100 (1 - x) (1 + y'[x]^2)^(3/2) + (1 + y'[x]^2)/y[x], 
  y[0] == 1, y[1] == 1}, y[x], {x, 0, 1}]

gives only the following message:

Power::infy: Infinite expression 1/0. encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`. >>

Is it because something is wrong with the ODE itself? or I will have to use some special options in NDSolve in order to handle such an ODE?

Answer 1

EDIT #2 (12.05.16)

I am greateful to the author of the OP for pointing out an error I made in EDIT #1 (which was also made by himself and user drN).

This leads me to revise my revision (EDIT #1) and return to the statement of the original solution: There is no solution to the original problem.

But can we can some words more on the question why there is no solution and "how far away from solubility" is the original system.

Experimenting as suggested in the first solution we find that the second derivative is so strong that it bends the curve up across the line y = 1 already for x < 1.

Hence we can possibly reach the goal y(x=1)=1 by reducing the impact of the second derivative by changing the formulas.

There are, of course, many ways to accomplish this.

We consider two possibilities

1) Reduce impact of second derivative by reducing the power of (1+y'[x]^2)

We already did the replacement 3/2 -> 1/2 "unnoticed" in EDIT #1. So these results can be reinterpreted in the current framework and are still of interest.

Letting 3/2->p we find for a critical power pc = 0.5917 defined such that there is no soluton for p > pc

With[{p = 0.5917}, 
  yy[x_] = y[x] /. 
    NDSolve[{y''[x] == 
        204/100 (1 - x) (1 + y'[x]^2)^p + (1 + y'[x]^2)/y[x], 
       y[0] == 1, y[1] == 1}, y[x], {x, 0, 1.5}][[1]]];
Plot[{1, yy[x]}, {x, 0, 1}] (* graph not shown here *)

{yy[0], yy[1]}
{1., 1.}

2) Reduce impact of second derivative by reducing the multiplicative factor "a" instead of 204/100

By trial and error we find a critical value ac = 1.3107 defined such that for a>ac there is no solution possible.

With[{a = 1.3107}, 
  yy[x_] = y[x] /. 
    NDSolve[{y''[x] ==(*204/100 *)
        1.3107 (1 - x) (1 + y'[x]^2)^(3/2) + (1 + y'[x]^2)/y[x], y[0] == 1, 
       y[1] == 1}, y[x], {x, 0, 1.5}][[1]]];
Plot[{1, yy[x]}, {x, 0, 1}] (* graph not shown here *)

{yy[0], yy[1]}
{1., 1.}

EDIT #1 (26.03.16)

Let us apply the shooting method, i.e. solve the equation with initial conditions and look for the value of y[1].

Let y[0] = 1 and y'[0] equal to some value v which we shall vary in order to hopefully get y[1] = 1:

Revising my previous statement we find that there are even two different solutons to the problem, i.e. the solution is not defined uniquely by the boundary conditions. This meakes the problem even more interesting.

I was mistaken in thinking of a monotoneous dependence of y[1] of y'[0] (=v). Closer inspection reveals that this was wrong.

Here's the code to experiment with:

v = -10; yy[x_] = 
 y[x] /. NDSolve[{y''[x] == 
      204/100 (1 - x) (1 + y'[x]^2)^(3/2) + (1 + y'[x]^2)/y[x], y[0] == 1, 
     y'[0] == v}, y[x], {x, 0, 2}][[1]];
Plot[{1, yy[x]}, {x, 0, 2}]

For a systematic approach let uns define the value y1 = y[1] as a function of the slope v at x = 0

y1[v_] := (y[x] /. 
    NDSolve[{y''[x] == 
        204/100  (1 -  x) (1 + y'[x]^2)^(1/2) + (1 + y'[x]^2)/y[x], y[0] == 1,
        y'[0] == v}, y[x], {x, 0, 2}][[1]]) /. x -> 1

The plot shows that the two values v ~= -3.5 and v ~= -2.6 make y[1] == 1:

Plot[y1[v] - 1, {v, -4, -2}, AxesLabel -> {"y'[0]", "y[1]-1"}]

Greater precision is obtained manually adjusting v such that y1 = 1.

For solution 1

v = -2.582496; 
yy1[x_] = y[x] /. 
  NDSolve[{y''[x] == 
      204/100  (1 -  x) (1 + y'[x]^2)^(1/2) + (1 + y'[x]^2)/y[x], y[0] == 1, 
     y'[0] == v}, y[x], {x, 0, 2}][[1]];
yy1[1]

(* Out[491]= 0.99999 *)

For solution 2:

v = -3.59259; 
yy2[x_] = y[x] /. 
  NDSolve[{y''[x] == 
      204/100  (1 -  x) (1 + y'[x]^2)^(1/2) + (1 + y'[x]^2)/y[x], y[0] == 1, 
     y'[0] == v}, y[x], {x, 0, 2}][[1]];
yy2[1]

(* Out[496]= 1. *)

Plotting both solutons together

Plot[{1, yy1[x], yy2[x]}, {x, 0, 1.2}, 
 PlotLabel -> "The two solutons of the boundary value problem", 
 AxesLabel -> {"x", "y[x]"}, PlotRange -> {0, 1.5}]

Original solution

Thus boundary value problem has no soluton.

Here we go to show it:

Let us solve the equation with initial conditions and look for the value of y1 obtained.

Let y[0] = 1 and y'[0] equal to some value v which we shall vary in order to hopefully get y1 == 1:

Here's the code to experiment with:

v = -10; yy[x_] = 
 y[x] /. NDSolve[{y''[x] == 
      204/100 (1 - x) (1 + y'[x]^2)^(3/2) + (1 + y'[x]^2)/y[x], y[0] == 1, 
     y'[0] == v}, y[x], {x, 0, 2}][[1]];
Plot[{1, yy[x]}, {x, 0, 2}]