# NDSolve unable to solve nonlinear differential equation

I have an action $$S$$,

$$S = \int dx \frac{1}{z^d} \sqrt{1 + \frac{z'(x)^2}{f(z)}}$$

where the Lagrangian $$L$$ is,

$$L = \frac{1}{z^d} \sqrt{1 + \frac{z'(x)^2}{f(z)}}$$

I want to plot $$z(x)$$ and solve the value of $$S$$, however, I encountered some error. The equations of motion (EOM) are nonlinear and I know they can be solved using some numerical approximation techniques, actually, this post is related to (275496) where @AlexTrounev used the wavelet approach. The Lagrangian there is much more complicated than what is here, however, I found that they produced similar errors using NDSolve.

Actually, the Lagrangian here has a conserved quantity which allowed me to simplify $$S$$ by hand and solve it with the help of Mathematica. The reason I'm posting this is that I'm hoping that NDSolve can fully solve the EOM without resorting to conservation quantities since the Lagrangian is much simpler but that didn't work. I want to know if it is an NDSolve issue where it cannot find the correct technique (thus needing the wavelet approach) or an equation issue.

I want to emphasize that I really want to know what part of the equation is causing the issue (which was not cleared in the other post) so that I can understand it better.

The plot of $$z(x)$$ should have the form shown in a sample below (disregarding the scale/values on the axes),

Needs["VariationalMethods"]
f = 1 - (z[x]/zh)^(d + 1);
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;
eqeuler = EulerEquations[L, z[x], x];

wp = 20;
equations = SetPrecision[{eqeuler /. {d -> 3, zh -> 10}, z[10^-5] == 10^-5, z'[10] == 0}, wp];
eq = equations[[1]];
bc1 = equations[[2]];
bc2 = equations[[3]];

sols = NDSolve[{eq, bc1, bc2}, z, {x, 10^-5, 10}, Method -> "StiffnessSwitching", WorkingPrecision -> wp]

Power::infy: Infinite expression 1/0^3 encountered.
Power::infy: Infinite expression 1/0^4 encountered.
Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered.
Power::infy: Infinite expression 1/0^3 encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.
Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered.
Infinity::indet: Indeterminate expression 0 ComplexInfinity ComplexInfinity encountered.
General::stop: Further output of Infinity::indet will be suppressed during this calculation.
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 1.20.*^-5.

Plot[Evaluate[z[x] /. sols], {x, 10^-5, 10}, AspectRatio -> 3/4, Frame -> True, FrameLabel -> {"x", "z"}, PlotStyle -> {Blue, 24}, LabelStyle -> {Black, 20}, PlotRange -> Full, ImageSize -> Large, AspectRatio -> 3/4]


Update I have tried the shooting method as advised in the posts (237117) (72725), however, I still encounter problems like step size is effectively zero; singularity or stiff system suspected. I'm not sure if shooting method is the best way.

• if you change z'[10] == 0 to z'[10^-5] == 0 so that now both initial conditions are at same point, then this error goes away. (but you still get NDSolve::ndsz: At x == 0.000014311849291849864, step size is effectively zero; singularity or stiff system suspected. this tells me the problem is because you have initial conditions at different points? Commented Nov 11, 2022 at 17:20
• @Nasser This is a boundary condition problem, I cannot just simply do that change, the form of $z(x)$ should look as in the image that I posted (see update). Commented Nov 11, 2022 at 17:27
• I did not say to change the BC, I was simply pointing out an observation that the error goes away when both conditions are at same point, that is all. So the problem seems to be due to having the two conditions at different points. Commented Nov 11, 2022 at 17:30
• @Nasser Yes, you're right. So does this tell us that NDSolve cannot solve it as is, and maybe some sort of boundary value problem techniques are required? Commented Nov 11, 2022 at 18:07
• You can try to implement a shooting method, as shown in my answer here: mathematica.stackexchange.com/a/274790/61267 Commented Nov 11, 2022 at 19:50

We can solve this problem with "ExplicitEuler" as follows. First, we define equations

Needs["VariationalMethods"];
f = 1 - (z[x]/zh)^(d + 1);
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;
eqeuler = EulerEquations[L, z[x], x];

eq = eqeuler /. {d -> 3, zh -> 10}; bc = {z[10^-5] == 10^-5,
v[10] == 0};

s = Solve[eq, z''[x]][[1]] // Simplify;
eq01 = {v'[x] == s[[1, 2]] /. z'[x] -> v[x], z'[x] == v[x]};



Second, we define function

zm[e_?NumericQ] :=
Module[{eps = e, s1},
s1 = NDSolveValue[{eq01, z[10] == eps, v[10] == 0},
z[10^-5], {x, 10^-5, 10}, Method -> "ExplicitEuler",
StartingStepSize -> 10^-4]; s1];


Third, we compute root

sol =
FindRoot[zm[e] == 10^-5, {e, 991/100}, WorkingPrecision -> 30] //
Quiet

(*Out[]= {e -> 9.90980851516520491227795116651}*)


Check that sol gives initial condition with error of 10^-10

zm[e] /. sol

Out[]= 0.0000100001


Finally, we compute and plot solution

Z = NDSolveValue[{eq01, z[10] == e /. sol, v[10] == 0},
z, {x, 10^-5, 10}, Method -> "ExplicitEuler",
StartingStepSize -> 10^-4];

Plot[Z[x], {x, 10^-5, 10}, Frame -> True, PlotRange -> All]


• When I removed Quiet in FindRoot it produced errors like General::ovfl: Overflow occurred in computation., NDSolveValue::nlnum, InterpolatingFunction::dmval, I want to understand why it was produced and how to resolve them. Why do you just use Quiet? Commented Nov 18, 2022 at 11:00
• This message means that during computation FindRoot has tested region $zm<10^{-5}$. It is why we need to check root zm[e] /. sol. Commented Nov 18, 2022 at 12:10
• Is there any way to avoid this issue or anything to do to remedy this issue of FindRoot`? Commented Nov 20, 2022 at 9:46
• Yes, we can avoid, but algorithm is not so simple like above. Commented Nov 21, 2022 at 3:31