shooting method and stiffness problem for NDSolve

Question

I'm trying to solve numerically the following differential equation:

1/r[x]^5 + k/Sqrt[1 + r'[x]^2] - (k r[x]r''[x])/(1 + r'[x]^2)^(3/2) == 0

I can set boundary condition in x=x0 and x=xF with x0 and xF being the boundaries of the domain of the x variable, i.e., x=[x0,xF]. The particular values are not important, I can set what I want as long as xF>x0>0. I do not have constraint on the value of the first derivative at the boundaries.

This equation has a singularity for r''[x] in x=0. To avoid this problem, I have to use shooting method with appropriate boundary conditions.

sf = NDSolve[{1/r[x]^5 + k/Sqrt[1 + r'[x]^2] - (k r[x]r''[x])/(1 + r'[x]^2)^(3/2) == 0, r[0] == 1, r[1] == 10}, r, {x, 0.1, 100}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {r[0] == 1, r'[0] == 0}}, 
  MaxSteps -> 500]

If I do so, I obtain the following errors:

NDSolve::mxst: Maximum number of 200 steps reached at the point x == 0.48535728275604967`. >>
NDSolve::mxst: Maximum number of 200 steps reached at the point x == 0.9112853584099116`. >>
NDSolve::mxst: Maximum number of 200 steps reached at the point x == 0.9766603163139572`. >>
General::stop: Further output of NDSolve::mxst will be suppressed during this calculation. >>
FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option. >>
NDSolve::berr: There are significant errors {0.0395308,-0.00050132} in the boundary value residuals. Returning the best solution found. >>

But anyway I obtain a solution sf in the form of an interpolating function:

{{r -> InterpolatingFunction[{{0.1, 0.920942}},<>]}}

It is this solution reliable?

I think that the the differential equation may be affected by stiffness issues, so I tried to add as a sub Method options to solve stiffness method as indicated in the help:

sf = NDSolve[{1/r[x]^5 + k/Sqrt[1 + r'[x]^2] - (k r[x]r''[x])/(1 + r'[x]^2)^(3/2) == 0, r[0] == 1, r[1] == 10}, r, {x, 0.1, 100}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {r[0] == 1, r'[0] == 0},Method -> {"StiffnessSwitching", "NonstiffTest" -> False}}, 
  MaxSteps -> 500]

The problem is that I obtain several errors, and a discontinuous solution that it's clearly wrong.

I will really appreciate any help.

Samir

Answer 1

I do get some warning with your code but not same as yours, maybe you misremembered the k you chose?

Since you mentioned you may choose various parameters for the equation in the future, I'd like to point out that, though I don't know the exact theoretical explanation(maybe something about chaos theory?), based on my experience, choosing a proper initial value is the most important thing and seems to be the only way to avoid failures when using shooting method, which can be a hard work sometimes, but not for the specific example you give: with a slight change of the "StartingInitialConditions" I got a result without any warning:

k = 1;
eqn = 1/r[x]^5 + k/Sqrt[1 + r'[x]^2] - (k r[x] r''[x])/(1 + r'[x]^2)^(3/2);
sf = NDSolve[{eqn == 0, r[0] == 1, r[1] == 10}, r, {x, 0, 100}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {r[0] == 1, r'[0] == 1}}, 
   WorkingPrecision -> 16];

Plot[r[x] /. sf, {x, 0, 100}, PlotRange -> All, PlotLabel -> "Result"]
Plot[eqn /. sf, {x, 0, 100}, PlotLabel -> "Error Check"]