3

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

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
sam84
  • 497
  • 6
  • 15
  • A first remark is about the k parameter you have in there. A numerical solution is not possible for parametric equations. – tchronis Dec 15 '13 at 16:50
  • I forget to specify. At present I fixed k to any numerical integer value greater that zero. In the above case k=1. In the future I will have to determine the value of the parameter k from a set of data points, for which the solution of the differential equation is the modelling. – sam84 Dec 16 '13 at 09:14

1 Answers1

4

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

enter image description here enter image description here

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • You are right, initial conditions are crucial and with r[0]==1, r'[0]==1, the solution is easily obtained. But if for instance I put the constraint r[10]==10, instead of r[1]=10, the solution is not able anymore to satisfy the constrain I impose. – sam84 Dec 18 '13 at 17:33
  • In fact the complete problem is 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[3] == 25, r[149] == 45}, r, {x, 3, 149}, Method -> {"Shooting", "StartingInitialConditions" -> {r[3] == 25, r'[3] == 1}}, WorkingPrecision -> 16]; But the solution gives me r[3] /. sf =105.2684511703500 – sam84 Dec 18 '13 at 17:35
  • 1
    @sam84 I believe what you need is still a better initial condition, though I can't find it now :). As I mentioned above, you may have a hard time searching a proper initial condition. See this and this for example. – xzczd Dec 19 '13 at 03:02