# The numerical solution of a nonlinear ODE: boundary value problem

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 == 1, y == 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?

• Could you share a little about the physics that this equation is solving? That way someone may be able to tell you about any well established boundary condition issue or specification with physical insight, perhaps? It would seem that it is the final term in your equation (1 + Derivative[y][x]^2)/y[x] causing the issue. Mar 26, 2016 at 12:49
• It's from a Chinese article on determining the curve shape of a one dimensional liquid surface with specific surface tension and pressure difference at the surface. So you will see something from the nonlinear ODE that is similar to and actually derived from the curvature radius of the liquid profile curve. Mar 28, 2016 at 3:11
• Could you please cite the paper by any chance? Thank you. Mar 28, 2016 at 12:50

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 == 1, y == 1}, y[x], {x, 0, 1.5}][]];
Plot[{1, yy[x]}, {x, 0, 1}] (* graph not shown here *)

{yy, yy}
{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 == 1,
y == 1}, y[x], {x, 0, 1.5}][]];
Plot[{1, yy[x]}, {x, 0, 1}] (* graph not shown here *)

{yy, yy}
{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.

Let y = 1 and y' equal to some value v which we shall vary in order to hopefully get y = 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 of y' (=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 == 1,
y' == v}, y[x], {x, 0, 2}][];
Plot[{1, yy[x]}, {x, 0, 2}] For a systematic approach let uns define the value y1 = y 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 == 1,
y' == v}, y[x], {x, 0, 2}][]) /. x -> 1


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

Plot[y1[v] - 1, {v, -4, -2}, AxesLabel -> {"y'", "y-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 == 1,
y' == v}, y[x], {x, 0, 2}][];
yy1

(* Out= 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 == 1,
y' == v}, y[x], {x, 0, 2}][];
yy2

(* Out= 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 = 1 and y' 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 == 1,
y' == v}, y[x], {x, 0, 2}][];
Plot[{1, yy[x]}, {x, 0, 2}] • Could this be a version related issue? I just started a fresh kernel and the problem solved fine "right out of the box" using NDSolveValue: NDSolveValue[{(y^\[Prime]\[Prime])[x] == 204/100 (1 - x) (1 + Derivative[y][x]^2)^(1/2 ) + ( 1 + Derivative[y][x]^2)/y[x], y == 1, y == 1}, y[x], {x, 0, 1}] Mar 26, 2016 at 13:09
• @drN You wrote 1/2 instead of 3/2. In this case it works out fine also in my version which is 10.1.0 for Microsoft Windows (64-bit) (March 24, 2015). But the problem was a different one. Mar 26, 2016 at 13:32
• Thank you! I found by using Newton's method it is possible to find the two initial conditions more accurate, with different start points: (1)FindRoot[obj[x] == 1, {x, -2}, Evaluated -> False] // NumberForm[#, 15] &  (1) FindRoot[obj[x] == 1, {x, -4}, Evaluated -> False] // NumberForm[#, 15] &`. The two start points were of course easily be identified by plotting the curve as in the answer. thanks Mar 28, 2016 at 3:22
• @user6043040 I am not aware of the literature on the topic, but I would start with en.wikipedia.org/wiki/Shooting_method Mar 28, 2016 at 8:03
• @user6043040 You are right, I changed the exponent to 1/2 because it was done in the comments of drN and you and I didn't check it. I shall make a brief EDIT in my answer to summaize the status. May 12, 2016 at 14:23