1
$\begingroup$

I am trying to evaluate:

FindRoot[Xfnew[γ, w, Hi, mx]*92*10^23*γ^(4/(1 + w)) == 1198/10000, {Hi, hs}, 
    MaxIterations -> 10000]

for some starting value hs, and γ and w are known numbers listed below in the definitions. I have to evaluate this FindRoot in a range of mx, and plot it, something like:

LogLogPlot[Hi /. FindRoot[Xfnew[γ, w, Hi, mx]*9.2*10^24*γ^(4/(1 + w)) == 
    0.1198, {Hi, 1.6946*10^-8}], {mx, 10^(-4), 10^(-10)}, PlotStyle -> Green]

I have successfully evaluated this in the range of mx (0.1 to 0.001), and (10^-10 to 10^-19), with starting values 0.0046 and 0.1 respectively.

But in the range 10^(-4) to 10^(-10), I am not able to. I tried using a "for loop" and break this range to form a list {mx,Hi} with small steps of mx(and doing a ListLogLogPlot afterwards), and perturbing the initial starting of FindRoot, like this:

For[i = 1, i < 2000, i++,h[0] = hs; h[i] = Hi /. Check[FindRoot[
    Xfnew[γ, w, Hi, mxs - i*step]*92.*10^23*γ^(4/(
    1 + w)) == 1198./10000, {Hi, h[i - 1]}, MaxIterations -> 10000], 
    {Hi -> 25}]; If [h[i] == 25, Break[], AppendTo[hin, {mxs - i*step, h[i]}]]]

where hs and mxs are some starting values of Hi and mx, and I break the loop if FindRoot generates error messages. And "step" is some stepsize I have chosen.

But, currently my step size for this loop has decreased so much, that it would take me ages to cover the entire range.I decrease my stepsize, whenever I get the error:

Machine precision is insufficient to achieve the accuracy

Is there a faster way to cover this range? I have tried increasing Working precision, it still doesnot help.

Function definitions

Xfnew is a result of a complicated list of functions listed below, and the values of gamma and w are listed below:: (Xfnew and the values of gamma and w are listed at the End) (I apologize for the Greek letters)

g = 61;
k2 = 25/100;    
k1 = 2/10;
ai = 1;
Trh[γ_?NumericQ, Hi_?NumericQ] := k2*γ*(( Hi)^(1/2))    
arh[γ_?NumericQ, w_?NumericQ] := ai*γ^(-4/(3 (1 + w)))
Hrh[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] := 
    Hi (a/ai)^(-3 (1 + w)/2)
H[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] := 
    Hi γ^2 (a/arh[γ, w])^-2
Trh11[a_?NumericQ, γ_?NumericQ, w_?NumericQ, 
    Hi_?NumericQ] := ((k1 (γ Hi)^(1/2))/((1 + 3/5 w)^((1/
    4)))) ((a^(-3 (1 - w)/2) - a^-4)^((1/4))) 
T[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] := 
    Trh[γ, Hi] arh[γ, w]/a
nXeq[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, 
    mx_?NumericQ, T_Symbol] := g/(2 π^2)
    mx^2 T[a, γ, w, Hi] BesselK[2, mx/T[a, γ, w, Hi]]
σv0[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
    mx_?NumericQ, T_Symbol] := 3 π ((mx^2)/
    8) (3 BesselK[1, mx/T[a, γ, w, Hi]]^2/(5 BesselK[2, mx/T[a, γ, w, Hi]]^2) + 2/5 + 
    4/5 T[a, γ, w, Hi]/mx BesselK[1, mx/T[a, γ, w, Hi]]/
    BesselK[2, mx/T[a, γ, w, Hi]] + 8/5 T[a, γ, w, Hi]^2/mx^2) 
σv1[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
    mx_?NumericQ, T_Symbol] := π (
    mx^2) (11/20 BesselK[1, mx/T[a, γ, w, Hi]]^2/
    BesselK[2, mx/T[a, γ, w, Hi]]^2 + 9/20 + 13/20 T[a, γ, w, Hi]/
    mx BesselK[1, mx/T[a, γ, w, Hi]]/BesselK[2, mx/T[a, γ, w, Hi]] + 
    13/10 T[a, γ, w, Hi]^2/mx^2)   
σvtotal[a_?NumericQ, γ_?NumericQ, w_?NumericQ, 
    Hi_?NumericQ, mx_?NumericQ, T_Symbol] := 
    4 σv0[a, γ, w, Hi, mx, T] + (45 + 12) σv1[a, γ, w, Hi, mx, T]
g1[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, 
    mx_?NumericQ] := (mx*a^2)/(Hrh[a, γ, w, Hi]*Trh[γ, Hi]^3) σvtotal[
    a, γ, w, Hi, mx, Trh11] nXeq[a, γ, w, Hi, mx, Trh11]^2; 
g2[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, 
    mx_?NumericQ] := (mx*a^2)/(H[a, γ, w, Hi]*Trh[γ, Hi]^3) σvtotal[
    a, γ, w, Hi, mx, T] nXeq[a, γ, w, Hi, mx, T]^2; 
Xfnew[γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, mx_?NumericQ] := 
    NIntegrate[g1[a, γ, w, Hi, mx], {a, 1, arh[γ, w]}] + 
    NIntegrate[g2[a, γ, w, Hi, mx], {a, arh[γ, w], ∞}]
γ = 1/10;
w = 5/10;
$\endgroup$
2
$\begingroup$

Experimentation indicates that a larger value of WorkingPrecision indeed is necessary. However, too large a value is computationally expensive and also sometimes causes the computation to crash. WorkingPrecision -> 30 seems to be a good compromise. In addition, computations for mx in the range of {10^-6, 10^-4} do not converge well without very good initial guesses. An initial guess as a function of mx fit to values of Hi for mx in the vicinity of 10^-2 and of 10^-7 and again a bit of experimentation suggests 0.2 mx^1.7 for mx > 10^-6. Then,

ParallelTable[{mx = N[10^i], 
    Hi /. FindRoot[Xfnew[γ, w, Hi, mx]*92*10^23*γ^(4/(1 + w)) == 1198/10000, 
    {Hi, Max[10^-11, .2 mx^(1.7)]}, WorkingPrecision -> 30, PrecisionGoal -> 6]}, 
    {i, -19, -1, .1}];
ListLogLogPlot[%, Joined -> True, PlotStyle -> Green, ImageSize -> Large, 
    AxesLabel -> Automatic, LabelStyle -> {Black, Bold, 14}]

enter image description here

$\endgroup$

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.