0
$\begingroup$

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.

Improve this question
$\endgroup$
1
$\begingroup$

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]

fig1

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. fig2

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$ – DiSp0sablE_H3r0 Jul 30 '18 at 9:33
  • $\begingroup$ 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. $\endgroup$ – Alex Trounev Jul 30 '18 at 11:38
  • $\begingroup$ Thanks for the suggestions and the added comments. Really helpful. Cheers!!! $\endgroup$ – DiSp0sablE_H3r0 Jul 30 '18 at 20:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.