3
$\begingroup$

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

s = NDSolve[{
      (z + 2 r[z] r'[z]) (1 + r'[z]^2) -z r[z] r''[z] == 0, 
      r[0.01] == 0.0001, 
      r'[0.01] == 10
     }, 
     r, {z, 0.01, 10}
    ]

But, when I set the initial condition to $z=0$, there are singularities, so I tried to set it as 0.01. Also, for $r(z=0.01)=0$ there is a singularity, so I set $r(z=0.01)=0.0001$.

However, when I try to solve the equation for any value that I put in the boundary conditions, I'm getting "step size is effectively zero; singularity or stiff system suspected." Is there any method to sort this out?

Improve this question
$\endgroup$
2
  • $\begingroup$ The singularity is clear: this equation algebraically implies $r'' = (1+r'^2)\frac{1}{r}$ plus other stuff not involving $r$. Because $1+r'^2\ge 1$ (assuming you want non-Complex solutions), it would be hard for $r$ to approach $0$ while its second derivative is blowing up like $1/r$. This indicates you will have extreme trouble finding any solution for small values of $r$. $\endgroup$
    – whuber
    Commented Oct 3, 2012 at 15:56
  • $\begingroup$ @William: Can you give more context of the problem for trying out altered boundary conditions? $\endgroup$
    – Narasimham
    Commented May 28, 2020 at 6:54

2 Answers 2

Reset to default
1
$\begingroup$

Problem seems to arise from initial values chosen.

Symbols $(r-z)$ used could represent a meridian of a surface of revolution comprising dimensionless terms.

$$ \dfrac {rr''}{1+r^{'2}}= 2 \dfrac{rr'}{z} + 1 $$

Qualitatively, this ODE is a combination of two tendencies. Separately taken,

the first ODE

$$ \dfrac {rr''}{1+r^{'2}}=k$$

represents curvature ratio progressively (1,2,3,4) ( catenary, cone, sphere and Weingarten profiles ). Integrating one obtains slope

$$ \cos \phi = {r}/{r_{(max/min)}} = k $$

and the second ODE

$$ \dfrac{rr'}{z} =\dfrac{rr'}{zz'} =const\; i.e.,\; \dfrac{dr^2}{dz^2} =\dfrac{2rr'}{2zz'}=const$$ can be put in a form to represent conicoids $ r^2 + k z^2 =1 $ that transition from hyperbolic paraboloid, cone, ellipsoid surfaces of revolution.

How they are transitioning be seen in this interesting problem, and it was a pleasure to realize some of its details.

Boundary conditions are chosen to capture full real solution for these sets.

The initial high slope and low radius tried out earlier seems not a suitable choice of initial conditions. In view of ellipsoidal closed profiling tendency expected from above logic I took zero starting initial value slope for one profiles and also varied start slope from a fixed initial point and another profiles set. Terminated the profile on crash by ode stiffness, infinite curvatures as happens for cusps. The surface is influenced only slightly due to action of $z$ term.

Profiles in Transition 

zmin = -4; zmax = 2.205;
NDSolve[{(1 + 2 r[z] r'[z]/z) - r[z] r''[z]/(1 + r'[z]^2) == 0, 
   r[zmin] == 1, r'[zmin] == 0}, r, {z, zmin, zmax}];
rn[u_] = r[u] /. First[%];
g3 = Plot[rn[z], {z, zmin, 2.205}, GridLines -> Automatic, 
   PlotStyle -> Thick];
CurvRatio3 = 
  Plot[rn[z] rn''[z]/(1 + rn'[z]^2), {z, zmin, 2.205}, 
   GridLines -> Automatic, PlotLabel -> Curvtr_Ratio, 
   PlotStyle -> Thick];
Plot[{rn[z], rn'[z], rn''[z]}, {z, zmin, 2.205}, 
 GridLines -> Automatic, PlotStyle -> Thick]
zmin = -2.5; zmax = 1.326;
NDSolve[{(1 + 2 r[z] r'[z]/z) - r[z] r''[z]/(1 + r'[z]^2) == 0, 
   r[zmin] == 1, r'[zmin] == 0}, r, {z, zmin, zmax}];
rn[u_] = r[u] /. First[%];
g4 = Plot[rn[z], {z, zmin, zmax}, GridLines -> Automatic, 
   PlotStyle -> Thick];
CurvRatio4 = 
  Plot[rn[z] rn''[z]/(1 + rn'[z]^2), {z, zmin, zmax}, 
   GridLines -> Automatic, PlotLabel -> Curvtr_Ratio, 
   PlotStyle -> Thick];
zmin = -1; zmax = 0.6731;
NDSolve[{(1 + 2 r[z] r'[z]/z) - r[z] r''[z]/(1 + r'[z]^2) == 0, 
   r[zmin] == 1, r'[zmin] == 0}, r, {z, zmin, zmax}];
rn[u_] = r[u] /. First[%];
g5 = Plot[rn[z], {z, zmin, 0.6731}, GridLines -> Automatic, 
   PlotStyle -> Thick];
CurvRatio5 = Plot[rn[z] rn''[z]/(1 + rn'[z]^2), {z, zmin, 0.6731}];
Show[{g3, g4, g5}, AspectRatio -> 0.5, PlotRange -> {{-4, 4}, {0, 4}}]
Show[{CurvRatio3, CurvRatio4, CurvRatio5}, AspectRatio -> .8, 
 PlotRange -> {{-4, 4}, {-2.5, 1.1}}]
ParametricPlot3D [{rn[z] Cos[t], z, rn[z] Sin[t]}, {z, zmin, 
  0.6731}, {t, -Pi/2,  Pi/2}, PlotStyle -> Yellow, Axes -> False , 
 Boxed -> False]

Appreciate your comments, especially about the unified but transitioning nature of profiles from one non-linear DE to the other in this manner.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for this. This is a problem I posted about back in 2012. It was related to some holographic calculations in the context of AdS/CFT. $\endgroup$
    – Wiliam
    Commented May 29, 2020 at 16:39
2
$\begingroup$

This problem seems to be really stiff as whuber has pointed out. Using the BDF method and allowing NDSolve to ruminate about the order, didn't quite change the solution (z value where stiffness was encountered).

NDSolve[{(z + 2 r[z] r'[z]) (1 + r'[z]^2) - z r[z] r''[z] == 0, 
  r[0.01] == 0.0001, r'[0.01] == 10}, r, {z, 0.01, 10}, 
 MaxStepFraction -> 1/101, 
 Method -> {"BDF", "MaxDifferenceOrder" -> 5}]

Try these commands out as pointed out by several authors here and here:

These answers helped me out with issues I had with stiff equations.

You can look around and find heaps more information.

?NDSolve`BDF
?NDSolve`LSODA
Improve this answer
$\endgroup$
2
  • $\begingroup$ Huh, I haven't said anything about these equations. You intended to address whuber? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 3, 2012 at 16:14
  • $\begingroup$ @J.M. Yes, I meant whuber... I think I got that wrong. $\endgroup$
    – dearN
    Commented Oct 3, 2012 at 20:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.