# Problem with a numeric 1/0 in NDSolveValue

The following code tries to solve three coupled differential equations with some fixed parameters

\$PrePrint = # /. {Csc[z_] :> 1/Defer@Sin[z],
Sec[z_] :> 1/Defer@Cos[z], Cot[z_] :> Defer@Cos[z]/Defer@Sin[z],
Csch[z_] :> 1/Defer@Sinh[z], Sech[z_] :> 1/Defer@Cosh[z],
Coth[z_] :> Defer@Cosh[z]/Defer@Sinh[z]} &;


These are the parameters

    ϵ := 1/100
L := 10
d := 1
ϕ := 0
χm := 5 π/12
χp := 4 π/12
αp := 1
N5 := 1
M5 := 0
N3 := 100
ΔN3 := 0
gYM := 1
δ =
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)])];

α = -(M5/4) Cosh[δ] +
Sqrt[(π^2 N3)/gYM^2 + M5^2/16 Cosh[δ]^2];

αh = (gYM^2 α)/(4 π);


These are some functions with the parameters inside. I took these parameters as much simple I can

h1 = αp (-I α Sinh[v] -
M5/4 Log[
Tanh[(I π)/4 - (v - δ)/
2]]) + αp (I α Sinh[vb] -
M5/4 Log[Tanh[-((I π)/4) - (vb - δ)/2]]) //
FullSimplify;
h2 = αp αh (Cosh[v] + Cosh[vb]) // FullSimplify;
w = D[D[h1 h2, vb], v] // FullSimplify;
F1 = 2 h1 h2 D[h1, v] D[h1, vb] - h1^2 w // FullSimplify;
F2 = 2 h1 h2 D[h2, v] D[h2, vb] - h2^2 w;
subv = {v -> x[x2] + I y[x2], vb -> x[x2] - I y[x2]};
f42 = 2 ((F1 F2)/w^2)^(1/4) /. subv // FullSimplify // PowerExpand
(* 20 Sqrt[π] Cosh[x[x2]]^2 *)

ρ2 =
4 (F1 F2 w^2)^(1/4)/(h1 h2) /. subv /.
Sin[2 y[x2]] :> 2 Cos[y[x2]] Sin[y[x2]] // FullSimplify //
PowerExpand
(* 20 Sqrt[π] *)


This is the lagrangian and the diffrential equations involved

Lag = Factor@
FullSimplify@f42 (u'[x2]^2/u[x2]^2 + 2/
u[x2]^2) + ρ2 (x'[x2]^2 + y'[x2]^2);

equ = Factor@FullSimplify@(D[Lag, u[x2]] - D[D[Lag, u'[x2]], x2]);

eqx = Factor@FullSimplify@(D[Lag, x[x2]] - D[D[Lag, x'[x2]], x2]);

eqy = Factor@FullSimplify@(D[Lag, y[x2]] - D[D[Lag, y'[x2]], x2]);

pdes = {equ == 0, eqx == 0, eqy == 0};


These are the boundary conditions

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;

Ld = L + d Sin[ϕ] // N;
L - d Sin[ϕ] > 0;

in = {u0, x0, y0} // N;
out = {u1, x1, y1} // N;

{x20, x21} // N;


Here I recollect the boundary conditions

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

bcsin = {u[x20] == u0, x[x20] == x0, y[x20] == y0};

bcsd = {bcsin, u' == 0, x' == 0, y' == 0} // Flatten;


Here I try to solve the equations

sols = NDSolve[{pdes, bcs}, {u, x, y}, {x2, x20, x21}];


Here I plot the solutions

Plot[sols[][x2], {x2, x20, x21}]
Plot[sols[][x2], {x2, x20, x21}]
Plot[sols[][x2], {x2, x20, x21}]


This code produces the following errors

Power::infy: Infinite expression 1/0. encountered. >>

Power::infy: Infinite expression 1/0.^2 encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Power::infy: Infinite expression 1/0. 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 encountered. >>

General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at x2 == -1.. >>

And, when I try to plot

NDSolve::dsvar: -0.999959 cannot be used as a variable. >>

NDSolve::dsvar: -0.959143 cannot be used as a variable. >>

NDSolve::dsvar: -0.918326 cannot be used as a variable. >>

General::stop: Further output of NDSolve::dsvar will be suppressed during this calculation. >>

Notice that even the equations are 2nd order, I give only fixed points (BVP). I think the shooting method cannot be applied when I have more than one equations.

• Must be that you divide by zero... somewhere in your convoluted code. – Henrik Schumacher Mar 2 '18 at 18:58
• @HenrikSchumacher I understand that. The problem is where it could happen, because the equations involved, if you see, are quite simple. – Patrick El Pollo Mar 2 '18 at 19:02
• @corey979 Those arent variables, those are the extreme values of x2, I mean x2=[x20,x21] – Patrick El Pollo Mar 2 '18 at 19:27

This BVP can be solved with using Haar wavelets. The system of equations and boundary conditions after simplification have a form

{pdes, bcs} = {{2 - Derivative[u][x2]^2 + 2*Tanh[x[x2]]*u[x2]*Derivative[u][x2]*Derivative[x][x2] + u[x2]*Derivative[u][x2] == 0,
-2*Sinh[2*x[x2]] - Sinh[2*x[x2]]*Derivative[u][x2]^2 + 2*u[x2]^2*Derivative[x][x2] == 0, Derivative[y][x2] == 0},
{u[-1] == Sqrt/100, x[-1] == ArcSinh, y[-1] == Pi/12, u == Sqrt/100, x == ArcSinh, y == Pi/6}};


We change variables as x2->-1+2 x2 and {f, q, t} -> {u, x, y}, so the new variable is in {x2,0,1}. Code with Haar wavelets for 16 colocation points is

L = 2; A = 0; B = 1; J = 3; M = 2^J; dx = (B - A)/(2*M); h1[x_] := Piecewise[{{1, A <= x <= B}, {0, True}}]; p1[x_, n_] := (1/n!)*(x - A)^n;
h[x_, k_, m_] := Piecewise[{{1, Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]}, {-1, Inequality[(1 + 2*k)/(2*m), LessEqual, x, Less, (1 + k)/m]}}, 0]
p[x_, k_, m_, n_] := Piecewise[{{0, x < k/m}, {(-(k/m) + x)^n/n!, Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]},
{((-(k/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, (1 + 2*k)/(2*m) <= x <= (1 + k)/m},
{((-(k/m) + x)^n + (-((1 + k)/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, x > (1 + k)/m}}, 0]
xl = Table[A + l*dx, {l, 0, 2*M}]; xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2*M + 1}];
f2[x_] := Sum[af[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*h1[x];
f1[x_] := Sum[af[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*p1[x, 1] + f10;
f0[x_] := Sum[af[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*p1[x, 2] + f10*x + f00;
q2[x_] := Sum[aq[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + e0*h1[x];
q1[x_] := Sum[aq[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + e0*p1[x, 1] + q10;
q0[x_] := Sum[aq[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + e0*p1[x, 2] + q10*x + q00;
t2[x_] := Sum[at[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + c0*h1[x];
t1[x_] := Sum[at[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + c0*p1[x, 1] + t10;
t0[x_] := Sum[at[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + c0*p1[x, 2] + t10*x + t00;
bc1 = {f0 == Sqrt/100, q0 == ArcSinh, t0 == Pi/12};
bc2 = {f0 == Sqrt/100, q0 == ArcSinh, t0 == Pi/6};
var = Flatten[Table[{af[i, j], at[i, j], aq[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
varM = Join[{a0, c0, e0, f10, f00, q10, q00, t10, t00}, var];
eqf[x2_] := 2*L^2 - f1[x2]^2 + 2*Tanh[q0[x2]]*f0[x2]*f1[x2]*q1[x2] + f0[x2]*f2[x2];
eqq[x2_] := -2*Sinh[2*q0[x2]]*L^2 - Sinh[2*q0[x2]]*f1[x2]^2 + 2*f0[x2]^2*q2[x2];
eqt[x2_] := t2[x2];
eq = Flatten[Table[{eqq[x] == 0, eqf[x] == 0, eqt[x] == 0}, {x, xcol}]];
eqM = Join[eq, bc1, bc2];
sol = FindRoot[eqM, Table[{varM[[i]], 0.1}, {i, Length[varM]}], MaxIterations -> 10000]; lstf = Table[{-1 + L*x, Evaluate[f0[x] /. sol]}, {x, 0, 1, 0.01}];
lstq = Table[{-1 + L*x, Evaluate[q0[x] /. sol]}, {x, 0, 1, 0.01}];
lstt = Table[{-1 + L*x, Evaluate[t0[x] /. sol]}, {x, 0, 1, 0.01}];


Numerical solution visualization

{ListLinePlot[{lstf, lstq, lstt}, PlotLegends -> {"u", "x", "y"}],
ListLinePlot[lstf, AxesLabel -> {"x2", "u"}]} 