Singularity and/or non-significant digits in BVP with NDSolve

I have the following code, in which I try to solve a system of three coupled differential equation numerically:

The code:

    (*Parameters*)

(*precision:20,50,\[Infinity]*)

pr := \[Infinity]
prr := pr - 5

cut := 1/1000000000
\[Epsilon] := SetPrecision[1/10, pr]
L := SetPrecision[1, pr]
d := SetPrecision[1/100, pr]
\[Phi] := SetPrecision[\[Pi]/8, pr]
\[Chi]m := SetPrecision[0, pr]
\[Chi]p := SetPrecision[0, pr]
\[Alpha]p := SetPrecision[1, pr]
N5 := SetPrecision[2, pr]
N3 := SetPrecision[20, pr]
\[CapitalDelta]N3 := SetPrecision[2, pr]
gYM := SetPrecision[10, pr]
\[Delta] =
SetPrecision[
1/2 Log[1/(
gYM^2 N5^2 (2 N3 - \[CapitalDelta]N3)) (2 gYM^2 N3 N5^2 +
4 \[Pi]^2 \[CapitalDelta]N3^2 +
Sqrt[(2 gYM^2 N3 N5^2 + 4 \[Pi]^2 \[CapitalDelta]N3^2)^2 -
gYM^4 N5^4 (4 N3^2 - \[CapitalDelta]N3^2)])], pr];
\[Alpha] =
SetPrecision[-(N5/4) Cosh[\[Delta]] +
Sqrt[(\[Pi]^2 N3)/gYM^2 + N5^2/16 Cosh[\[Delta]]^2], pr];
\[Alpha]h = SetPrecision[(gYM^2 \[Alpha])/(4 \[Pi]), pr];

(*harmonic functions and definitions*)

h1 = \[Alpha]p (-I \[Alpha] Sinh[v] -
N5/4 Log[
Tanh[(I \[Pi])/4 - (v - \[Delta])/
2]]) + \[Alpha]p (I \[Alpha] Sinh[vb] -
N5/4 Log[Tanh[-((I \[Pi])/4) - (vb - \[Delta])/2]]);
h2 = \[Alpha]p \[Alpha]h (Cosh[v] + Cosh[vb]);
w = D[D[h1 h2, vb], v];
F1 = 2 h1 h2 D[h1, v] D[h1, vb] - h1^2 w;
F2 = 2 h1 h2 D[h2, v] D[h2, vb] - h2^2 w;
f42s = N[2 ((F1 F2)/w^2)^(1/4), prr];
\[Rho]2s = N[2 h2^2 ((F1 w^2)/F2^3)^(1/4), prr];
subv = {v -> x[x2] + I y[x2], vb -> x[x2] - I y[x2]};

(*Coeff functions inside the eoms*)

(*THESE FUNCTIONS ARE INSIDE THE DIFFERENTIAL EQUATIONS BELOW*)

(*these functions will be evaluated at the initial values given \
below, BUT they give undefined values or non-significant digits \
depending on the precision*)

Logfx = SetPrecision[D[Log[f42s /. subv], x[x2]], prr];
Logfy = SetPrecision[D[Log[f42s /. subv], y[x2]], prr];
Log\[Rho]x = SetPrecision[D[Log[\[Rho]2s /. subv], x[x2]], prr];
Log\[Rho]y = SetPrecision[D[Log[\[Rho]2s /. subv], y[x2]], prr];
f\[Rho]x =
SetPrecision[D[f42s /. subv, x[x2]]/(\[Rho]2s /. subv), prr];
f\[Rho]y =
SetPrecision[D[f42s /. subv, y[x2]]/(\[Rho]2s /. subv), prr];

(*Lagrangian and equations of motion*)

(*f42 is Subscript[f, 4](x,y)^2 and \[Rho]2 is \[Rho](x,y)^2*)
(*Subscript[f, 4](x,y) and \[Rho](x,y) are given*)

(*SUBSTITUTIONS*)

subeq = { f42[x[x2], y[x2]] -> f42,
\!$$\*SuperscriptBox[\(f42$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x[x2], y[x2]] -> f42y,
\!$$\*SuperscriptBox[\(f42$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x[x2], y[x2]] ->
f42x, \[Rho]2[x[x2], y[x2]] -> \[Rho]2,
\!$$\*SuperscriptBox[\(\[Rho]2$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x[x2], y[x2]] -> \[Rho]2x,
\!$$\*SuperscriptBox[\(\[Rho]2$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x[x2], y[x2]] -> \[Rho]2y};

(*I dont know if the last line will be clear as code, so I added a \
copy of it below*)

(*subeq={ f42[x[x2],y[x2]]\[Rule]f42,D[ \
f42[x[x2],y[x2]],y[x2]]\[Rule]f42y,D[ \
f42[x[x2],y[x2]],x[x2]]\[Rule]f42x,\[Rho]2[x[x2],y[x2]]\[Rule]\      [Rho]2,\
D[\[Rho]2[x[x2],y[x2]],x[x2]]\[Rule]\[Rho]2x,D[\[Rho]2[x[x2],y[x2]],y[\
x2]]\[Rule]\[Rho]2y};*)

subu = {f42x -> A f42, f42y -> B f42};

subxy = {f42x -> F \[Rho]2,
f42y -> G \[Rho]2, \[Rho]2x -> H \[Rho]2, \[Rho]2y -> J \[Rho]2};

subwarp2 = {A -> Logfx, B -> Logfy, H -> Log\[Rho]x, J -> Log\[Rho]y,
F -> f\[Rho]x, G -> f\[Rho]y};

(*Lagrangian*)

Lag = f42[x[x2], y[x2]] (u'[x2]^2/u[x2]^2 + 2/u[x2]^2) + \[Rho]2[
x[x2], y[x2]] (x'[x2]^2 + y'[x2]^2);

(*SET OF EQUATIONS OF MOTION where A,B,F,G,H and J are defined above*)

(*eom for u(x2)*)

equ =  (u[x2]^2)/(2 f42) (D[Lag, u[x2]] - D[D[Lag, u'[x2]], x2]) /.
subeq /. subu // Expand;

(*eom for x(x2)*)

eqx =  1/(2 \[Rho]2) (D[Lag, x[x2]] - D[D[Lag, x'[x2]], x2]) /.
subeq /. subxy // Expand;

(*eom for y(x2)*)

eqy =  1/(2 \[Rho]2) (D[Lag, y[x2]] - D[D[Lag, y'[x2]], x2]) /.
subeq /. subxy // Expand;

(*collecting the equations*)

equu = (equ + u''[x2]) /. subwarp2;
eqxx = (eqx + x''[x2]) /. subwarp2;
eqyy = (eqy + y''[x2]) /. subwarp2;
pdes = {u''[x2] == equu, x''[x2] == eqxx, y''[x2] == eqyy};

(*Boundary conditions*)

(*initial and final values of x2*)

x20 = SetPrecision[-d Cos[\[Phi]], prr];
x21 = SetPrecision[d Cos[\[Phi]], prr];

(*initial values of u,x and y*)

u0 = SetPrecision[\[Epsilon] Sqrt[
1 + ((L - d Sin[\[Phi]])/\[Epsilon])^2], prr];
x0 = SetPrecision[ArcSinh[(L - d Sin[\[Phi]])/\[Epsilon]], prr];
y0 = SetPrecision[\[Pi]/2 - \[Chi]m, prr];

(*Final values*)

u1 = SetPrecision[\[Epsilon] Sqrt[
1 + ((L + d Sin[\[Phi]])/\[Epsilon])^2], prr];
x1 = SetPrecision[ArcSinh[(L + d Sin[\[Phi]])/\[Epsilon]], prr];
y1 = SetPrecision[\[Pi]/2 - \[Chi]p, prr];
(*required condition*)
L - d Sin[\[Phi]] > 0

(*boundary conditions*)

bcs = SetPrecision[{x[x20] == x0, u[x20] == u0, y[x20] == y0,
x[x21] == x1, u[x21] == u1, y[x21] == y1}, prr];

(*slight deformation of the initial conditions*)

cor := SetPrecision[1/1000000, prr]
bcs2 = SetPrecision[{x[x20] == x0 + cor, u[x20] == u0 + cor,
y[x20] == y0 + cor, x[x21] == x1 + cor, u[x21] == u1 + cor,
y[x21] == y1 + cor}, prr];

(*Solving eoms*)

(*NDSolve[Flatten[{pdes,bcs}],{u[x2],x[x2],y[x2]},{x2,x20,x21},Method\
\[Rule]{"ExplicitRungeKutta","DifferenceOrder"\[Rule]2},\
WorkingPrecision\[Rule]70][[1]];*)
subin = SetPrecision[{u[x20] -> u0, x[x20] -> x0, y[x20] -> y0}, prr];
subin2 = SetPrecision[{u[x20] -> u0 + cor, x[x20] -> x0 + cor,
y[x20] -> y0 + cor}, prr];

(*Evaluating the eoms at the initial value of x20*)

bcsin = {x[x20] == x0, u[x20] == u0, y[x20] == y0};
bcsin2 = {x[x20] == x0 + cor, u[x20] == u0 + cor, y[x20] == y0 + cor};
(*EVALUATION AT THE ORIGINAL INITIAL VALUES*)

NDSolve[Flatten[{pdes, bcs}], {u[x2], x[x2], y[x2]}, {x2, x20, x21},
Method -> {"Shooting", "StartingInitialConditions" -> bcsin},
WorkingPrecision -> MachinePrecision][[1]];
NDSolve::ndsz: At x2 == -0.0092388, step size is effectively zero; singularity or stiff system suspected. >>

(*EVALUATION AT THE "DEFORMED INITIAL VALUES*)

NDSolve[Flatten[{pdes, bcs2}], {u[x2], x[x2], y[x2]}, {x2, x20, x21},
Method -> {"Shooting", "StartingInitialConditions" -> bcsin2},
WorkingPrecision -> MachinePrecision][[1]];

NDSolve::ndsz: At x2 == -0.0092388, step size is effectively zero; singularity or stiff system suspected. >>


Now the problem. Initially, I set the variable to some numerical values, and there is a chain of definitions involving them, ok. Since some of them are trigonometric functions, square roots, etc; I set infinite precision (which displays numbers inside these functions, e.g Sin[pi/2] instead Sin[3.14.../2]), so coefficient functions are finally evaluated in those parameters and depend only on functions and numbers (with infinite precision, i.e non-functional expressions).

I know that Mathematica deal BVP with Shooting method. I also know that NDSolve needs numerical values inside the equations, that is why I set WorkingPrecision to MachinePrecision. The problema appear when I run the code and get

NDSolve::ndsz: "At x0 == (number), step size is effectively zero; singularity or stiff system suspected."


So the problems looks to be at the beginning. I evaluated the coefficient functions inside the PDEs and saw that indeed x=x0 gives something like

"Infinite expression 1/(0.*10^-37+0.*10^-37\I) encountered."


Maybe there are analytic zeroes inside the coefficient functions, I thought, so I applied Chop to set zero small numbers. When I evaluate the coefficient functions I get

no significant digits to display


And something like 0.x10^-20 or 0.x10^50. These "numbers" spoil the evaluation and the method, I think... producing something like 1/0 in the middle steps.

Any suggestion? I hope you understand the problem.

• Please display the actual code that produces the error. Otherwise, it is nearly impossible for the reader to provide a useful answer. Commented Nov 2, 2017 at 2:07
• @bbgodfrey see below. Commented Nov 2, 2017 at 22:44
• Please edit your question to include the code in your "Answer", which surely will be deleted, because it is not, in fact, an answer. When you do so, replace **Parameters** etc. by (*Parameters*). Commented Nov 2, 2017 at 23:00
• @bbgodfrey Done Commented Nov 2, 2017 at 23:20
• I have eliminated the error messages listed in the question but need more information to provide good guesses for the Shooting method. What do you expect the solution to look like? Is it real everywhere? Is it approximately symmetric about the origin? Are derivatives of the functions approximately zero at x2 == 0? Etc. Commented Nov 3, 2017 at 2:06

The error message given in the question occurs because the ODEs are singular at y[x2] == Pi/2, which is the boundary condition given for y (see discussion at the end of this answer). SetPrecision has no effect on this problem and only clutters the code. With it and other unnecessary material deleted, the code becomes

ϵ = 1/10; L = 1; d = 1/100; ϕ = π/8; χm = 1/1000; χp = 1/1000; αp = 1; N5 = 2;
N3 = 20; ΔN3 = 2; gYM = 10;
δ = 1/2 Log[1/(gYM^2 N5^2 (2 N3 - ΔN3)) (2 gYM^2 N3 N5^2 + 4 π^2 ΔN3^2 +
Sqrt[(2 gYM^2 N3 N5^2 + 4 π^2 ΔN3^2)^2 - gYM^4 N5^4 (4 N3^2 - ΔN3^2)])];
α = -(N5/4) Cosh[δ] + Sqrt[(π^2 N3)/gYM^2 + N5^2/16 Cosh[δ]^2];
αh = (gYM^2 α)/(4 π);

h1 = αp (-I α Sinh[v] - N5/4 Log[Tanh[(I π)/4 - (v - δ)/2]]) +
αp (I α Sinh[vb] - N5/4 Log[Tanh[-((I π)/4) - (vb - δ)/2]]);
h2 = αp αh (Cosh[v] + Cosh[vb]);
w = D[D[h1 h2, vb], v];
F1 = 2 h1 h2 D[h1, v] D[h1, vb] - h1^2 w;
F2 = 2 h1 h2 D[h2, v] D[h2, vb] - h2^2 w;
f42s = 2 ((F1 F2)/w^2)^(1/4);
ρ2s = 2 h2^2 ((F1 w^2)/F2^3)^(1/4);
subv = {v -> x[x2] + I y[x2], vb -> x[x2] - I y[x2]};

Logfx = D[Log[f42s /. subv], x[x2]];
Logfy = D[Log[f42s /. subv], y[x2]];
Logρx = D[Log[ρ2s /. subv], x[x2]];
Logρy = D[Log[ρ2s /. subv], y[x2]];
fρx = D[f42s /. subv, x[x2]]/(ρ2s /. subv);
fρy = D[f42s /. subv, y[x2]]/(ρ2s /. subv);

subeq = {f42[x[x2], y[x2]] -> f42, D[f42[x[x2], y[x2]], y[x2]] -> f42y,
D[f42[x[x2], y[x2]], x[x2]] -> f42x, ρ2[x[x2], y[x2]] -> ρ2,
D[ρ2[x[x2], y[x2]], x[x2]] -> ρ2x, D[ρ2[x[x2], y[x2]], y[x2]] -> ρ2y};
subu = {f42x -> A f42, f42y -> B f42};
subxy = {f42x -> F ρ2, f42y -> G ρ2, ρ2x -> H ρ2, ρ2y -> J ρ2};
subwarp2 = {A -> Logfx, B -> Logfy, H -> Logρx, J -> Logρy, F -> fρx, G -> fρy};

Lag = f42[x[x2], y[x2]] (u'[x2]^2/u[x2]^2 + 2/u[x2]^2) +
ρ2[x[x2], y[x2]] (x'[x2]^2 + y'[x2]^2);
equ = (u[x2]^2)/(2 f42) (D[Lag, u[x2]] - D[D[Lag, u'[x2]], x2]) /.
subeq /. subu // Expand;
eqx =  1/(2 ρ2) (D[Lag, x[x2]] - D[D[Lag, x'[x2]], x2]) /. subeq /. subxy // Expand;
eqy =  1/(2 ρ2) (D[Lag, y[x2]] - D[D[Lag, y'[x2]], x2]) /. subeq /. subxy // Expand;
pdes = {equ == 0, eqx == 0, eqy == 0} /. subwarp2;

x20 = -d Cos[ϕ];
x21 = d Cos[ϕ];

u0 = ϵ Sqrt[1 + ((L - d Sin[ϕ])/ϵ)^2];
x0 = ArcSinh[(L - d Sin[ϕ])/ϵ];
y0 = π/2 - χm;

u1 = ϵ Sqrt[1 + ((L + d Sin[ϕ])/ϵ)^2];
x1 = ArcSinh[(L + d Sin[ϕ])/ϵ];
y1 = π/2 - χp;

bcs = {x[x20] == x0, u[x20] == u0, y[x20] == y0,
x[x21] == x1, u[x21] == u1, y[x21] == y1};

NDSolveValue[{pdes, bcs}, {u[x2], x[x2], y[x2]}, {x2, x20, x21}];


Unfortunately, Method -> "Shooting", which is automatic here, invariably failed, even with initial guesses explicitly given, due to failure to converge or, more commonly, division by zero. Note that χm = 1/1000; χp = 1/1000 are used in the code above to avoid singularities at the boundaries that occurred in the question. Some insight can be obtained from treating the calculation as an IVP instead of a BVP,

bcst = {u[x20] == u0, x[x20] == x0, y[x20] == y0, u'[x20] == 0,
x'[x20] == -10, y'[x20] == 10};
NDSolveValue[{pdes, bcst}, {u[x2], x[x2], y[x2]}, {x2, x20, x21}];
Plot[{Chop[%[[1]]], Chop[%[[2]]], Chop[%[[3]]]}, {x2, x20, x21}, ImageSize -> Large,
LabelStyle -> Directive[Black, Bold, Medium], AxesLabel -> {"x2", "u, x, y"}]


Large changes to x'[x20] and y'[x20] have almost no impact on the curves shown.

Comments by the OP indicate that the desired solution is expected to be approximately symmetric about x2 == 0. This suggests trying something like

bcst = {x[0] == x0 - .19, u[0] == u0, y[0] == Pi/2 - .08, x'[0] == 0,
y'[0] == 0, u'[0] == 0};
NDSolveValue[{pdes, bcst}, {u[x2], x[x2], y[x2]}, {x2, x20, x21}];
Plot[%, {x2, x20, x21}, ImageSize -> Large,
LabelStyle -> Directive[Black, Bold, Medium], AxesLabel -> {"x2", "u, x, y"}]
(%% /. x2 -> x20) // Chop


(* {1.00111, 2.98889, 1.36707} *)


It is easy in this way to approximately match the boundary conditions for x and u given in the question. However, the boundary conditions for y given in the question are far from being matched. Based on my many attempts to obtain even approximately the desired y boundary conditions, I believe that they cannot be achieved. (Remember that the ODEs are highly nonlinear.)

Discussion of singularity at y == Pi/2

The overall structure of the ODEs can be obtained as follows.

subwarphf = {A -> HoldForm@Logfx, B -> HoldForm@Logfy, H -> HoldForm@Logρx,
J -> HoldForm@Logρy, F -> HoldForm@fρx, G -> HoldForm@fρy}


However, the six quantities Logfx etc. have LeafCounts in the tens of thousands, and Mathematica cannot simplify them or expand them about y[x2] == Pi/2 in a reasonable amount of time. They can, however be plotted in minutes, showing that Logρy and fρx are singular there. From the structure of the eqx and eqy, it is evident that they too must be singular unless x' and y' vanish identically there.

Addendum: Solution Satisfying u and x Boundary Conditions

In comments below, the OP requested a solution with u and x satisfying the boundary conditions given in the question, and y satisfying 0 < y[x2] < Pi/2 but otherwise unconstrained. One such solution is given by

bcst = {u[x20] == u0, x[x20] == x0, y[x20] == 1.36707,
u'[x20] == 0.3504, x'[x20] == -18.9337, y'[x20] == 20.0155};
NDSolveValue[{pdes, bcst}, {u[x2], x[x2], y[x2]}, {x2, x20, x21}];
Plot[%, {x2, x20, x21}, ImageSize -> Large,
LabelStyle -> Directive[Black, Bold, Medium], AxesLabel -> {"x2", "u, x, y"}]
Plot[First@%%, {x2, x20, x21}, ImageSize -> Large,
LabelStyle -> Directive[Black, Bold, Medium], AxesLabel -> {"x2", "u"}]
{%%% /. x2 -> x20, %%% /. x2 -> x21} // Chop


(* {{1.00118, 2.99441, 1.36707}, {1.0088, 3.00202, 1.3649}} *)


Comparing u[x21] and x[x21] with

u1 // N
(* 1.0088 *)
x1 // N
(* 3.00202 *)


shows that the solution obtained here satisfies the x21 boundary conditions to five significant figures. I obtained the solution parameters starting from those for the second plot in the answer and then adjusting u'[x20] and x'[x20] to satisfy the two x21 boundary conditions. Using the "Shooting" Method, on the other hand, invariably failed. It seems likely, nonetheless, that using ParametricNDSolve combined with FindRoot, appropriately tuned would have yielded the answer I obtained by trial and error.

• Thanks for your answer. Indeed, y==Pi/2 as an initial and final value is not correct. Even taking Chi_m=Chi_p=1/1000 as you considered looks to give strange plots, sorry for not being clear in that detail. I'm trying differente values for y at x20 and x21, but this is matter of choosing a suitable set of numbers. Moreover, I dont understand the second set of conditions bcst which give the second plot, e.g u[0] == u0. Commented Nov 5, 2017 at 8:46
• bcst sets first derivatives of y, x, and y to zero at x2 == 0 and integrates from there to x20 and x21, thereby giving a solution which is symmetric about x0 = 0, as you suggested in a comment. I chose values for the three functions themselves at x2 == 0 to give values at the end points that were reasonably close to those in the question, except of course for y. If you planned to try different values for y at the endpoints, what were your criteria for a successful choice? Commented Nov 5, 2017 at 13:57
• The initial and endpoints of x y and u are x20 and x21, which depend on L, d and phi (which must combine into a positive exression givenin the code); also Chi_m, Chi_p and epsilon. My mistake was to consider the initial and final values of y as equal. A successful choice of parameters should be epsilon small, e.g. 1/10 or 1/100. About y, it must be greater than zero and less than Pi/2, and different at x20 and x21 (different Chis). It is good to choose vanishing derivative of x,y,u at x20==0. Commented Nov 5, 2017 at 20:32
• @resanrom If I understand you correctly, you wish u and x to have the endpoint values specified in the question, Pi/2 > y > 0 to have unequal values at the endpoints, and L - d Sin[ϕ] > 0. Is that correct? If so, then it is not possible to have first derivatives of u, x, and y all equal to zero at x2 == 0, although they could be close to zero there. Please clarify. Also, are there any other constraints on the desired solution? Commented Nov 5, 2017 at 21:17
• Yes, fixed points are specified before and unequal. Sorry again, I'm not being precise, I expect that x y and u behave as increasing functions, in general, until they reach x2==0, and then they decrease to their fixed value at x2==x21. In particular, I expect that u do that (like an "arc"), while for x and y I actually dont, but I could approximate that the also have vanishing first derivatives. No more constraints. Commented Nov 5, 2017 at 23:32