I have the following piece of code

(*UV cuttoff*)
rmax = 10;
(*Fixed input values*)
nc = 4;
nf = 4;
(*Value of α1*)
α1 = 0.95;

b0 = 1/(6 π) (11 nc - 2 (nc + 2) nf) ;

b1 =  1/(24 π^2) (34 nc^2 - 10 nc (nc - 2) nf  -  
     6 (nc + 1) (nc - 2) (nc - 2) nf) ;

eqalp1 = alp1'[mu1] + (b0 alp1[mu1]^2 + b1 alp1[mu1]^3) == 0;

sol1 = NDSolve [{ eqalp1, alp1[0] == α1}, 
   alp1, {mu1, 0, 10^24}];

alph1[mu1_] := First[Evaluate[alp1[mu1] /. sol1]]

α[mu_] := alph1[mu]

It works fine if I choose nf=1,2,3 but when I set nf=4 and higher I get the following error message


I tried to solve the differential equation not from zero but from other points. Every time I get a similar error. I also tried changing the tuned input for α1 in the beginning of the code, but the same message persists.

Can someone suggest something to bypass this?

P.S: I have used different combinations of max steps and working precision.I forgot to mention that previously. Finally, I tried to adopt this answer and use the WhenEvent command but it has no effect whatsoever.


Your ODE can be solved exactly, by the way. The time to reach a singularity is given by the following integral, if it converges:

 -1/(b0 alp1[mu1]^2 + b1 alp1[mu1]^3), 
 {alp1[mu1], α1, Infinity}]

(*  0.451388  *)

Therefore the IVP has a singularity around mu1 = 0.451. This agrees with the numerical integration performed by NDSolve (for nf = 4).

(The integral comes from separating the variables. Given $${d\alpha\over d\mu} = v(\alpha)\,,$$ we get, starting from the initial condition $\alpha(0)=\alpha_1$, $$\int_{\alpha_1}^{\infty} {d\alpha \over v(\alpha)} = \int_0^M d\mu = M$$ The value of $M$, if the integral on the left converges, is the time $\mu = M$ at which $\alpha \rightarrow \infty$.)

  • $\begingroup$ Thanks for the reply. I need the numerical interpolation for practical reasons thought as I have to use it as input for other equations. $\endgroup$ Nov 3 '19 at 19:45
  • 1
    $\begingroup$ @A_user_with_NoName You're welcome. I guess I was thinking you could use the exact solution to verify the results returned by NDSolve. $\endgroup$
    – Michael E2
    Nov 3 '19 at 20:06
  • $\begingroup$ this makes sense, of course. I guess that since I missed something so obvious it is a good time to stop for today and start fresh tomorrow. Thanks again for your reply! $\endgroup$ Nov 3 '19 at 20:13

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.