# Second Order Non Linear Differential Equation

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?

• 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$. – whuber Oct 3 '12 at 15:56
• @William: Can you give more context of the problem for trying out altered boundary conditions? – Narasimham May 28 at 6:54

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} = const ;\; or\quad \dfrac{dr^2}{dz^2} =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. 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]


• 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. – Wiliam May 29 at 16:39

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.

?NDSolveBDF