How to test if a solution of NDSolve is positive

Question

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.

Any help/comments are welcome.

Answer 1

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.