# How to test if a solution of NDSolve is positive

I want to solve a differential equation numerically and determine one of the parameters in the D.E. The goal is to determine $M$ such that the function goes asymptotically to zero or to put it in a better way close to zero, without ever getting negative values.

What I want to solve is the following

$$f''(r) + \frac{6}{r} f'(r) - \frac{1}{r^2} \Big(\frac{1}{2}\Big)^2 f(r) + \frac{6}{r^2}f(r) + \frac{1}{r^2} \frac{1}{2} f(r) + \frac{M^2}{(r^2 + lf(r)^2)^2}f(r) = 0$$

In order to calculate the $lf(r)$ in the above, some calculations are needed. Below, I give the code and I am breaking it in two pieces. The first is what is needed to calculate $lf(r)$ and the second is my attempt to solve the D.E.

The first piece of the code. I don't have any questions here and this is just to evaluate $lf(r)$ in the above D.E.

rmax := 10^4;
nc := 2;
nf1 := 2;
nf2 := 4;
\[Alpha]1 := 0.7;
b0 := 1/(6 \[Pi]) (11 (nc + 1) - 2 nf1 - 2 nf2 (nc - 1));
b1 := 1/(24 \[Pi]^2) (34 (nc + 1)^2 - (10 (nc + 1) +
6*(2 nc + 1)/4) 1/2 2 nf1 - (10 (nc + 1) + 6 nc) (nc - 1) nf2);
b0;
b1;
eqalp1 = alp1'[mu1] + (b0 alp1[mu1]^2 + b1 alp1[mu1]^3) == 0;
eqalp1;
sol1 = NDSolve [{ eqalp1, alp1[0] == \[Alpha]1},
alp1, {mu1, 0, 10^24}];
alph1[mu1_] := First[Evaluate[alp1[mu1] /. sol1]]
alph1[mu1];
\[Alpha][mu_] := alph1[mu]
\[Gamma][mu1_] := 3/(2 \[Pi]) \[Alpha][mu1] ((2 nc + 1)/4 + nc)
\[Delta][mu1_] := -2 \[Gamma][mu1]
muv = 0;
ode := D[r^3 L'[r], r] - \[Delta][Log[Sqrt[r^2 + L[r]^2] ]] r L[r] == 0
solm[v_] := NDSolve[{ode, L[v ] == v,
L'[v] == 0}, L, {r, v, rmax}, Method -> {"StiffnessSwitching"},
MaxSteps -> Infinity, PrecisionGoal -> 13, AccuracyGoal -> 13] //
First
LinfBm[ v_?NumericQ ] := (L[rmax] /. solm[v]) - muv
mir = m /. FindRoot[LinfBm[m], {m, 1}]
lf[\[Rho]_] := Evaluate[L[\[Rho]] /. solm[ mir]]


And the second part, the one that performs the NDSolve.

diracpde =
D[f[r], {r, 2}] + 6/r D[f[r], r] - 1/r^2 (1/2)^2 f[r] +
6 /r^2 f[r] + 1/r^2 (1/2) f[r] + M^2/(r^2 + lf[r]^2)^2 f[r];
diracforndsolve = {diracpde == 0, f[mir] == 1, f'[mir] == 0};
sol1 = NDSolve[diracforndsolve /. M :> 2, f[r], {r, mir, rmax}];
Plot[f[r] /. sol1, {r, mir, rmax},
PlotRange -> {{mir, rmax}, {-1, 1}},
BaseStyle -> {18, FontFamily -> "Times New Roman"},
AxesLabel -> {"\[Rho]", "top(\[Rho])"}, PlotStyle -> {Thick, Red}]


As mentioned above, the goal is to change the value of $M$ and find the solution that goes to zero, but never crosses to negative values.

One way of doing so, but I don't really like it is to plot from zero to one instead of -1 to 1.

My question is, if there is a better way to test it.

I have tried the following

h[r_] := f[r] /. sol1
FullSimplify[ForAll[r, mir <= r <= rmax, h[r] > 0]]

h[r_] := f[r] /. sol1
Reduce[ForAll[r, mir <= r <= rmax, h[r] > 0]]

FullSimplify[ForAll[r, mir <= r <= rmax \[Implies] h[r] > 0]]


motivated by some answers I found in other threads, but they don't do the job.

Another question I have been having is the following. If I set a lower $rmax$, i.e if I set it to $10$, the D.E has a solution with the feature tht I want and can be easily found, but for the case above $10^4$ it seems that it does not, or that I have done something wrong.

We make a substitution:

In[1]:= f[r] = Exp[q[r]]; eq =
D[f[r], {r, 2}] + 6/r D[f[r], r] - 1/r^2 (1/2)^2 f[r] + 6/r^2 f[r] +
1/r^2 (1/2) f[r] + M^2/(r^2 + lf[r]^2)^2 f[r];
diracpde = eq*Exp[-q[r]] // FullSimplify

Out[2]= 25/(4 r^2) + M^2/(r^2 + lf[r]^2)^2 + (
6 Derivative[1][q][r])/r +
Derivative[1][q][r]^2 + (q^\[Prime]\[Prime])[r]


Then if the function $q(r)$ is real it will always be $f(r)>0$. Let us investigate the behavior of $q(r)$ for $r$ approaching infinity:

DSolve[(diracpde /. {M -> 0}) == 0, q[r], r]


We have a simple analytic solution:

{{q[r] -> C[2] - (5 Log[r])/2}}


Consequently, $q(x)$ tends to minus infinity at r -> Infinity, respectively, $f(x)$ tends to zero, remaining positive everywhere. This solution does not depend on M. Let us show how the solutions for different $M$ vary as a function of $r$.

rmax = 10^4;
nc = 2;
nf1 = 2;
nf2 = 4;
\[Alpha]1 = 0.7;
b0 = 1/(6 \[Pi]) (11 (nc + 1) - 2 nf1 - 2 nf2 (nc - 1));
b1 = 1/(24 \[Pi]^2) (34 (nc + 1)^2 - (10 (nc + 1) + 6*(2 nc + 1)/4) 1/
2 2 nf1 - (10 (nc + 1) + 6 nc) (nc - 1) nf2);
eqalp1 = alp1'[mu1] + (b0 alp1[mu1]^2 + b1 alp1[mu1]^3) == 0;
eqalp1;
sol1 = NDSolve[{eqalp1, alp1[0] == \[Alpha]1}, alp1, {mu1, 0, 10^24}];
alph1[mu1_] := First[Evaluate[alp1[mu1] /. sol1]]
\[Alpha][mu_] := alph1[mu]
\[Gamma][mu1_] := 3/(2 \[Pi]) \[Alpha][mu1] ((2 nc + 1)/4 + nc)
\[Delta][mu1_] := -2 \[Gamma][mu1]
muv = 0;
ode := D[r^3 L'[r], r] - \[Delta][Log[Sqrt[r^2 + L[r]^2]]] r L[r] == 0
solm[v_] :=
NDSolve[{ode, L[v] == v, L'[v] == 0}, L, {r, v, rmax},
Method -> {"StiffnessSwitching"}, MaxSteps -> Infinity,
PrecisionGoal -> 13, AccuracyGoal -> 13] // First
LinfBm[v_?NumericQ] := (L[rmax] /. solm[v]) - muv
mir = m /. FindRoot[LinfBm[m], {m, 1}]
lf[\[Rho]_] := Evaluate[L[\[Rho]] /. solm[mir]]
diracpde =
D[f[r], {r, 2}] + 6/r D[f[r], r] - 1/r^2 (1/2)^2 f[r] +
6/r^2 f[r] + 1/r^2 (1/2) f[r] + M^2/(r^2 + lf[r]^2)^2 f[r];
diracforndsolve = {diracpde == 0, f[mir] == 1, f'[mir] == 0};
solm = ParametricNDSolveValue[diracforndsolve, f, {r, mir, rmax}, {M}];
LogLogPlot[Table[Abs[solm[k][r]], {k, 0, 4, 1}], {r, mir, rmax},
PlotRange -> All, BaseStyle -> {18, FontFamily -> "Times New Roman"},
AxesLabel -> {"\[Rho]", "top(\[Rho])"}, PlotStyle -> {Thick, Red},
PlotPoints -> 100]


And so, we see here one continuous solution, which corresponds to $M=0$ and four discontinuous solutions, which correspond to $M=1,2,3,4$. These data illustrate our assertion, which is not yet a theorem, but which must be proved. If we integrate not $f(r)$, but $q(r)$, then the situation changes (these are all the foci of numerical integration).

f[r] = Exp[q[r]]; eq =
D[f[r], {r, 2}] + 6/r D[f[r], r] - 1/r^2 (1/2)^2 f[r] + 6/r^2 f[r] +
1/r^2 (1/2) f[r] + M^2/(r^2 + lf[r]^2)^2 f[r];
diracpde = eq*Exp[-q[r]] // FullSimplify;
diracforndsolve = {diracpde == 0, q[mir] == 0, q'[mir] == 0};
solm = ParametricNDSolveValue[diracforndsolve,
q, {r, mir, rmax}, {M}];
LogLogPlot[Table[-solm[M][r], {M, 0, 4, .5}],{r, mir, rmax},
PlotRange -> {10^-10, 10^10}, AxesLabel -> {"r", ""},
PlotStyle -> Orange, PlotPoints -> 100, PlotLabel -> "-q(r)"]


Of the 9 values of {M,0,4,.5}, five solutions $q(r)$ are discontinuous, and four are continuous; there exist continuous solutions for the interval $0\le M\le Mmax$. But this still needs to be proved.

• Thanks for your time, however, the DSolve only works for the case M=0. In that case of course the solution is independent of M, since M is vanishing. If I set M=1,2,... it does not work. As I said in the question the goal is to find M such that the solution is a curve that asymptotically goes to zero for some value of r, but never crosses to negative values, and this is why I need to find a way to test if the numerical solution is positive. Cheers – Konstantinos Jul 30 '18 at 9:33
• We showed that the asymptotics of $q(r)$ for $r$ approaching infinity does not depend on M. Unfortunately, for finite $r$ the solutions $q(r)$ are discontinuous, therefore one can not guarantee that $f(r)$ differs from zero. Thus, it is required to find $M$ such that the function $q(r)$ has no discontinuities for finite $r$. We indicated one solution for M = 0. We must show that there are no other solutions. Numerical methods are useless here. – Alex Trounev Jul 30 '18 at 11:38
• Thanks for the suggestions and the added comments. Really helpful. Cheers!!! – Konstantinos Jul 30 '18 at 20:50