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

Question

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.

Answer 1

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.