ndsz : step size is effectively zero; singularity or stiff system suspected

Question

This is the first time I ask a question. I have seen many solutions to this error and tried but they are not working.

Here is the code.

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["FunctionApproximations`"];

Clear["Global`*"];

e = 1.60217662*10^-19;
c = 299792458;
m = 9.10938356*10^-31;
ϵ0 = 8.854187817*10^-12;
A = 10^17;
h = 6.62607004*10^-34; hbar = 1.0545718*10^-34;
k = 1/ϵ0 (2 e/h^3) (4 Pi/3);
x0 = (3/(4 *Pi)*A*1.60439572*10^-24/(4*10^14))^(1/3)/100;

xf = 10^0;

s = NDSolve[{1/x D[x y[x], {x, 2}] == 
    k ((e y[x]/c)^2 + 2 m*e y[x])^1.5, y[x0] == 10^6, y[xf] == 10^-5},
   y[x], {x, x0, xf}]

xf should be infinity but NDSolve doesn't allow infinity so I put a relatively large number.

y[xf] is also needed to be zero at infinity. I put a relatively small number since I expected that y[x] is close to zero when x is approaching infinity.

Here is the error:

At x == 0.112455460174063`, step size is effectively zero; singularity or stiff system suspected. >>

Answer 1

Analysis of the equation

It often is convenient to non-dimensionalize ODEs arising in physics, and this case is no exception. Rescale x and y by {x -> xn xs, y[x] -> yn[xn] ys}. Apply these transformations to the ODE, bringing the scale factors from the left side of the ODE, (xs^2/ys), to the right.

eqlast = (xs^2/ys) k ((e y[x]/c)^2 + 2 m*e y[x])^(3/2) /. 
    y[x] -> yn[t] ys /. x -> xn xs /. t -> xn
(* (k xs^2 (2 e m ys yn[xn] + (e^2 ys^2 yn[xn]^2)/c^2)^(3/2))/ys *)

Now, choose the scale factors, xs and ys, to eliminate the multipliers of yn[xn] and yn[xn]^2.

eq1 = eqlast[[1 ;; 3]] (eqlast[[4, 1, 1]])^(3/2) /. yn[xn] -> 1 /. xn -> 1
(* (2 Sqrt[2] k xs^2 (e m ys)^(3/2))/ys *)
eq2 = eqlast[[1 ;; 3]] (eqlast[[4, 1, 2]])^(3/2) /. yn[xn] -> 1 /. xn -> 1
(* (k xs^2 ((e^2 ys^2)/c^2)^(3/2))/ys *)
norms = FullSimplify[Solve[eq1 == 1 && eq2 == 1 && k > 0 && e > 0 && m > 0 && c > 0, 
    {xs, ys}, Reals], k > 0 && e > 0 && m > 0 && c > 0] // Last
(* {xs -> Sqrt[1/(c e k)]/(2 m), ys -> (2 c^2 m)/e} *)

Inserting these expressions into the transformed ODE yields,

eq = 1/xn D[xn yn[xn], {xn, 2}] == Simplify[eqlast /. %, k > 0 && e > 0 && m > 0 && c > 0]
(* (2 yn'[xn] + xn yn''[xn])/xn == (yn[xn] (1 + yn[xn]))^(3/2) *)

which is somewhat easier to work with than the ODE in its original form. The question seeks a solution with very small y at large x. To obtain such an asymptotic solution, assume that y^2 is much smaller than y (in absolute value). Then, it is apparent that a solution exists equal to cf xn^a, a and cf constants.

Unevaluated[1/xn D[xn yn[xn], {xn, 2}] - (yn[xn] (1 + yn[xn]))^(3/2)] /. 
    {(1 + yn[xn]) -> 1, yn[xn] -> cf xn^a}
(* a (1 + a) cf xn^(-2 + a) - (cf xn^a)^(3/2) *)

For this expression to vanish, xn^(-2 + a) == xn^(3 a/2) clearly is required. In other words, a == -4. Then, cf is determined from

Simplify[% /. a -> -4, cf > 0 && xn > 0]
(* -(((-12 + Sqrt[cf]) cf)/xn^6) *)

So, cf == 144, and this particular solution is yn[xn] == 144/xn^4. As we shall see, this is not the most general solution at large xn, but it does appear to be a separatrix.

Numerical Solution

With the constants given in the question, the scale factors become,

norms
(* {xs -> 3.46944*10^-12, ys -> 1.022*10^6} *)

and the starting point of the integration, x0, normalized to xs, becomes (rationalized for future use in NDSolve)

x0n = Rationalize[x0/xs /. norms, 0]
(* 723020393/5483051 *)

or about 131.865, and the specified value of y[x0], normalized to ys, becomes

y0n = Rationalize[10^6/ys /. norms, 0]
(* 33547031/34284995 *)

or about 0.978476. Conveniently, x0n is well away from the singularity at xn == 0. The problem, now, is to choose yn'[x0n] to that yn[xn] approaches the separatrix, obtained above, at large xn. Employing a variant on the method employed to solve question 147207 yields the desired value of yn'[x0n].

xfn = 2000;
s = ParametricNDSolve[{eq, yn[x0n] == y0n, yn'[x0n] == yp0, 
    WhenEvent[Re@yn[xn] < 120/xfn^4 || yn'[xn] > 1/100, "StopIntegration"]}, 
    yn, {xn, x0n, xfn}, {yp0, wp}, WorkingPrecision -> wp, AccuracyGoal -> Infinity, 
    MaxSteps -> 20000];

WhenEvent[ ... ] is used here to terminate integration, whenever a numerical solution veers away from the separatrix. Doing so is desirable, because one hundred or more iterations at high WorkingPrecision typically are needed to obtain a desired solution. A small amount of experimentation is sufficient to determine that the initial value for yn'[x0n] is near -1.39. A better approximation then is obtained from

plt = Plot[Quiet@(yn[yp0, 15]["Domain"] /. s)[[1, 2]], {yp0, -1.306, -1.304}, 
    PlotRange -> All]

The improved estimate, given by the spike in the plot, is located within the range,

Cases[plt, Line[a__] :> a, Infinity] // Last;
Position[%, Max[Last /@ %]][[1, 1]];
est = SetPrecision[{First[%%[[% - 1]]], First[%%[[% + 1]]]}, 45]
(* {-1.30493508434085847547123648837441578507423401, 
    -1.30493385315771526222761167446151375770568848} *)

This already narrow range is reduced much further with

dom[yp0_?NumericQ] := Quiet@(yn[yp0, 45]["Domain"] /. s)[[1, 2]] 
gg[bl0_, bu0_] := Module[{bl = N[bl0, 45], bu = N[bu0, 45], bt, db}, 
    Do[bt = (bl + bu)/2; db = (bu - bl)/100; 
    If[dom[bt] - dom[bt - db] > 0, bl = bt, bu = bt], {i, 120}]; bt];

(Note that FixedPoint could be used here instead of Do but appears to offer no advantage, at least for the present computation.)

r2000 = gg @@ est
(yn[r2000, 45][xfn] /. s) - 120/xfn^4
(* -1.30493476864289837272478679314595905244023206 *)
(* 8.3663144286145423790813208786836901834*10^-17 *)

About 80 sec was used by my PC to obtain this result. Such very high precision usually is necessary when seeking a solution asymptotically approaching a separatrix. Below is a plot of the numerical result, in blue, and the separatrix, in orange.

LogPlot[{yn[r2000, 45][xn] /. s, 144/xn^4}, {xn, x0n, xfn}, 
    AxesLabel -> {"x/xs", "y/ys"}, LabelStyle -> Directive[Bold, 11], PlotRange -> All]

As expected (and desired), the numerical solution is quite near the separatrix at large xn, even when viewed on a log plot.

Answer 2

Note

Because of an arguable backslide around v11, the solution in this answer doesn't work since then. See comments below for more information.

According to my personal experience, when approximating boundary condition (b.c.) at infinity with a b.c. at a large distance, usually the large distance should not be too large, and when solving nonlinear boundary value problem, manual setting a good enough initial guess for shooting method is often necessary. Your problem fits well in these rules:

e = 1.60217662*10^-19;
c = 299792458;
m = 9.10938356*10^-31;
ϵ0 = 8.854187817*10^-12;
A = 10^17;
h = 6.62607004*10^-34; hbar = 1.0545718*10^-34;
k = 1/ϵ0 (2 e/h^3) (4 Pi/3);
x0 = (3/(4*Pi)*A*1.60439572*10^-24/(4*10^14))^(1/3)/100;

xf = 50 10^-10;

sol = NDSolve[{1/x D[x y[x], {x, 2}] == k ((e y[x]/c)^2 + 2 m*e y[x])^1.5, y[x0] == 10^6, 
   y[xf] == 0}, y, {x, x0, xf}, 
  Method -> {"Shooting", "StartingInitialConditions" -> {y[xf] == 0, y'[xf] == 0}}]

Plot[y[x] /. sol, {x, x0, xf}, PlotRange -> All]

---

The following is a comparision between the solution given by bbgodfrey and me. As one can see, the 2 solutions are fairly consistent:

bbgodfreyplot = LogPlot[{yn[r2000, 45][xn] /. s, 144/xn^4} // Evaluate, {xn, x0n, xfn}, 
       AxesLabel -> {"x/xs", "y/ys"}, LabelStyle -> Directive[Bold, 11], 
   PlotRange -> All];
  
xzczdplot = 
 Unevaluated@
   LogPlot[y[x xs]/ys /. sol, {x, x0/xs, xf/xs}, PlotRange -> All, 
    PlotStyle -> {Red, Dashed, Thick}] /. norms;

Show[xzczdplot, bbgodfreyplot]