1
$\begingroup$

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),

Image

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.

$\endgroup$
8
  • 1
    $\begingroup$ 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? $\endgroup$
    – Nasser
    Nov 11, 2022 at 17:20
  • $\begingroup$ @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). $\endgroup$
    – mathemania
    Nov 11, 2022 at 17:27
  • 1
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Nov 11, 2022 at 17:30
  • $\begingroup$ @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? $\endgroup$
    – mathemania
    Nov 11, 2022 at 18:07
  • $\begingroup$ You can try to implement a shooting method, as shown in my answer here: mathematica.stackexchange.com/a/274790/61267 $\endgroup$
    – Hans Olo
    Nov 11, 2022 at 19:50

1 Answer 1

1
$\begingroup$

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]

Figure 1

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – mathemania
    Nov 18, 2022 at 11:00
  • 1
    $\begingroup$ This message means that during computation FindRoot has tested region $zm<10^{-5}$. It is why we need to check root zm[e] /. sol. $\endgroup$ Nov 18, 2022 at 12:10
  • $\begingroup$ Is there any way to avoid this issue or anything to do to remedy this issue of FindRoot? $\endgroup$
    – mathemania
    Nov 20, 2022 at 9:46
  • $\begingroup$ Yes, we can avoid, but algorithm is not so simple like above. $\endgroup$ Nov 21, 2022 at 3:31

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.