Getting FindRoot::cvmit and NDSolve::berr when evaluating an NDSolve expression

Question

I'm want to use NDSolve to compare non-linear equations and linear aproximations of them.

I have run into trouble with this equation:

EI = 1;
w = 24;
L = 1;
paso = 0.1;

sol2 =  
  NDSolve[
    {(15*Derivative[1][y][x]^2*Derivative[2][y][x]^3)/
        (1 + Derivative[1][y][x]^2)^(7/2) - 
     (3*Derivative[2][y][x]^3)/(1 + Derivative[1][y][x]^2)^(5/2) - 
     (9*Derivative[1][y][x]*Derivative[2][y][x]*Derivative[3][y][x])/
        (1 + Derivative[1][y][x]^2)^(5/2) + 
     Derivative[4][y][x]/(1 + Derivative[1][y][x]^2)^(3/2) == w/EI, 
     y[0] == 0, 
     Derivative[1][y][0] == 0, 
     Derivative[2][y][L] == 0, 
     Derivative[3][y][L] == 0},
    y, {x, 0, 1}, MaxSteps -> 1000, 
    PrecisionGoal -> 30, 
    StartingStepSize -> paso, 
    Method -> {"FixedStep", Method -> "ExplicitEuler"}]

When I run the program, Mathematica give me these errors:

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.
NDSolve::berr: The scaled boundary value residual error of 1.7005377392707127`*^8 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

And when I plot the results I obtain this:

instead something like this

Can someone tell me how to fix the errors I am getting?

Answer 1

Original Answer

Due to the NDSolve options chosen in the question,

`StartingStepSize -> 0.1, Method -> {"FixedStep", Method -> "ExplicitEuler"}]

the results obtained there are not at all accurate. Better results can be obtained by using default options, but only for smaller w/EI, which is represented as c below for convenience. Specifically, the code used in this answer is

c = 5.999; L = 1;
eq = 15*y'[x]^2*y''[x]^3/(1 + y'[x]^2)^(7/2) - 3*y''[x]^3/(1 + y'[x]^2)^(5/2) - 
      9*y'[x]*y''[x]*y'''[x]/(1 + y'[x]^2)^(5/2) + y''''[x]/(1 + y'[x]^2)^(3/2) 
sol2 =  Flatten@NDSolve[{eq == c, y[0] == 0, y'[0] == 0, y''[L] == 0, y'''[L] == 0},
    {y[x], y'[x], y''[x], y'''[x]}, {x, 0, L}, Method -> {"Shooting", 
    "StartingInitialConditions" -> {y''[0] == c/2, y'''[0] == -c}, "MaxIterations" -> 500}]

{y[x], y'[x], y''[x], y'''[x]} /. sol2 /. x -> 0
{y[x], y'[x], y''[x], y'''[x]} /. sol2 /. x -> L
(* {-4.81756*10^-9, 2.74752*10^-8, 2.9995, -5.999} *)
(* {6.80004, 54.7706, 4.45309*10^-19, -2.22686*10^-19} *)

GraphicsGrid[{{Plot[y[x] /. sol2, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y}],
     Plot[y'[x] /. sol2, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y'}]},
    {Plot[y''[x] /. sol2, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y''}],
     Plot[y'''[x] /. sol2, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y'''}]}}, 
    ImageSize -> 600]

I have been unable to obtain solutions for c >= 6, even though

"StartingInitialConditions" -> {y''[0] == c/2, y'''[0] == -c}

seems to be a very good initial guess. Specifically, in attempting to obtain a solution for c == 6, I tried "MaxIterations" -> 1000 together with WorkingPrecision -> 30 without success. Note that performing the calculation with c == 5.99 (instead of 5.999) yields

{y[x], y'[x], y''[x], y'''[x]} /. sol2 /. x -> L
(* {4.10476, 17.299, 7.70213*10^-23, 5.26227*10^-28} *)

In other words, y[L] and y'[L] increase rapidly as c == 6 is approached. It may be that no solution exists at this value of c and beyond.

Reduced Order ODE

Because y[x] does not appear in the ODE, it can be reduced to third order by the substitution

eq3 = eq /. Derivative[n_][y] -> Derivative[n - 1][z]
(* 15*z[x]^2*z'[x]^3/(1 + z[x]^2)^(7/2) - 3*z'[x]^3/(1 + z[x]^2)^(5/2) - 
      9*z[x]*z'[x]*z''[x]/(1 + z[x]^2)^(5/2) + z'''[x]/(1 + z[x]^2)^(3/2) *)

which can be solved quickly for c much closer to 6, after which y[x] can be obtained from y'[x] == z[x].

c = 5.999999;
sol3 =  Flatten@NDSolve[{eq3 == c, z[0] == 0, z'[L] == 0, z''L] == 0}, 
    {z[x], z'[x], z''[x]}, {x, 0, L}, Method -> {"Shooting", 
    "StartingInitialConditions" -> {z'[0] == c/2, z''[0] == -c}}];
sy = Flatten@NDSolve[{y'[x] == (z[x] /. sol3), y[0] == 0}, y[x], {x, 0, L}];

Flatten[{y[x] /. sy, {z[x], z'[x], z''[x]} /. sol3} /. x -> 0]
Flatten[{y[x] /. sy, {z[x], z'[x], z''[x]} /. sol3} /. x -> L]
(* {0., 6.47988*10^-8, 3., -6.} *)
(* {27.1176, 2159.61, -1.26109*10^-26, 8.09076*10^-27} *)

GraphicsGrid[{{Plot[y[x] /. sy, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y}],
     Plot[z[x] /. sol3, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y'}]},
    {Plot[z'[x] /. sol3, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y''}],
     Plot[z''[x] /. sol3, {x, 0, L}, PlotRange -> All, AxesLabel -> {x, y'''}]}},
    ImageSize -> 600]

This result reinforces the earlier conclusion that solutions to the nonlinear ODE exist over {x, 0, 1} only for c < 6. To determine whether the solutions obtained are unique, it is convenient to search for other solutions by means of

sol4 =  Flatten@NDSolve[{eq3 == c, z[0] == 0, z'[L] == 0, z''[L] == 0}, 
    {z[x], z'[x], z''[x]}, {x, 0, L}, Method -> {"Shooting", 
    "StartingInitialConditions" -> {z[L] == 3500}, "MaxIterations" -> 500}];

because "StartingInitialConditions" is needed for only one boundary condition at x == L. I have found that any initial guess between 0 and 3500 yields the solution given above, and "Shooting" does not converge otherwise. Other values of c yield analogous results.

Note the scaling y and x by a (a constant factor), and c by a^3 leaves the ODE unchanged but with the upper limit integration set to a^(-1/3). This suggests that integrating beyond N[4^(-1/3)] (i.e., 0.629961) is not possible for c = 24, consistent with the comment above by Bill.

Finally, it should be noted that the third-order ODE is autonomous and, therefore, can be further reduced to a second-order ODE by the procedure described in the answer to 140833. The advantage of doing so seems small in the present case, however.