2
$\begingroup$

I have the following Mathematica Code:

c = 2.99792458*^8;
A = 4.38*^-11;
t0 = 8.57697*10^17;
mB = 2.91*^-27;
ργ = 4.64*^-31;
σT = 6.65*^-29;
eta0 = 1.82*^18;

FreeElectronFractionData = {{3000, 1.0829044`}, {2909.991`, 
    1.0827902`}, {2819.9821`, 1.0825628`}, {2729.9731`, 
    1.0818148`}, {2639.9641`, 1.0792593`}, {2549.9551`, 
    1.0709482`}, {2459.9462`, 1.0503868`}, {2369.9372`, 
    1.0224601`}, {2279.9282`, 1.0056397`}, {2189.9193`, 
    1.0009301`}, {2099.9103`, 1.0001211`}, {2009.9013`, 
    1.0000111`}, {1919.8923`, 0.9999931`}, {1829.8834`, 
    0.99996694`}, {1739.8744`, 0.99984236`}, {1649.8654`, 
    0.99910856`}, {1559.8565`, 0.99387158`}, {1469.8475`, 
    0.95045371`}, {1379.8385`, 0.69166608`}, {1289.8296`, 
    0.2628594`}, {1199.8206`, 0.059562751`}, {1109.8116`, 
    0.0097843423`}, {1019.8026`, 0.0013015781`}, {929.79367`, 
    0.00017836243`}, {839.7847`, 0.000036054186`}, {749.77573`, 
    0.000012683805`}, {659.76676`, 7.0873662`*^-6}, {569.75779`, 
    5.0157033`*^-6}, {479.74881`, 3.98173`*^-6}, {389.73984`, 
    3.377567`*^-6}, {299.73087`, 2.9942744`*^-6}, {209.7219`, 
    2.740874`*^-6}, {119.71293`, 2.5702004`*^-6}, {29.70396`, 
    2.4566939`*^-6}};

etaT[t_] := t*(2*c + A*t)/c
tZ[z_] := Sqrt[t0^2/(z + 1)]
etaZ[z_] := etaT[tZ[z]]
zEta[η_] := t0^2/((-c + Sqrt[c]*Sqrt[c + A*η])/A)^2 - 1
a[z_] := 1/(z + 1)

R[z_, ρB_] := (3*ρB)/((a[z]^3*(4*ργ))/a[z]^4)
freeElectronFraction = Interpolation[FreeElectronFractionData];
electronDensity[
  z_, ρB_] := (freeElectronFraction[z]*ρB*(1 + z)^3)/mB
scatterRate[η_, ρB_] := 
 electronDensity[zEta[η], ρB]*σT*a[zEta[η]]*c
τ[η_, ρB_] := 
 NIntegrate[scatterRate[etaPrime, ρB], {etaPrime, η, eta0}]
etaStar[ρB_] := η /. 
  Quiet[FindRoot[τ[η, ρB] == 1, {η, etaZ[1100]}]]
sStar[z_, ρB_] := 
  Integrate[
   c/Sqrt[3*(1 + R[zEta[etaPrime], ρB])], {etaPrime, 0, etaZ[z]}];
DLSS[z_] := (t0*(A*t0*z + 2*c*(1 + z - Sqrt[1 + z])))/(2 + z)
AngularScale[ρB_] := (z = zEta[etaStar[ρB]];
  sStar[z, ρB]/DLSS[z])
AngularScale[5.453*10^-27]

Plot[AngularScale[ρ], {ρ, 5*10^-27, 6*10^-27}]

NMinimize[{Abs[
   AngularScale[ρ] - 1.0411], ρ > 5*10^-27 && ρ < 
    6*10^-27}, {ρ}]

During evaluation of In[91]:= NIntegrate::nlim: etaPrime = \[Eta] is not a valid limit of integration.

During evaluation of In[91]:= NIntegrate::inumr: The integrand (1.36452*10^37 \[Rho] <<1>>[-1+1.41129*10^15/(<<1>>)^2])/(-2.99792*10^8+17314.5 Sqrt[2.99792*10^8+4.38*10^-11 etaPrime])^4 has evaluated to non-numerical values for all sampling points in the region with boundaries {{5.17952*10^16,5.37952*10^16}}.

During evaluation of In[91]:= NIntegrate::inumr: The integrand (1.36452*10^37 \[Rho] <<1>>[-1+1.41129*10^15/(<<1>>)^2])/(-2.99792*10^8+17314.5 Sqrt[2.99792*10^8+4.38*10^-11 etaPrime])^4 has evaluated to non-numerical values for all sampling points in the region with boundaries {{5.17952*10^16,5.37952*10^16}}.

During evaluation of In[91]:= NIntegrate::inumr: The integrand (1.36452*10^37 \[Rho] <<1>>[-1+1.41129*10^15/(<<1>>)^2])/(-2.99792*10^8+17314.5 Sqrt[2.99792*10^8+4.38*10^-11 etaPrime])^4 has evaluated to non-numerical values for all sampling points in the region with boundaries {{5.17952*10^16,5.37952*10^16}}.

During evaluation of In[91]:= General::stop: Further output of NIntegrate::inumr will be suppressed during this calculation.

During evaluation of In[91]:= FindRoot::nlnum: The function value {-1.+NIntegrate[scatterRate[etaPrime,\[Rho]],{etaPrime,\[Eta],eta0}]} is not a list of numbers with dimensions {1} at {\[Eta]} = {5.17952*10^16}.

During evaluation of In[91]:= ReplaceAll::reps: {FindRoot[\[Tau][\[Eta],\[Rho]]==1,{\[Eta],etaZ[1100]}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

During evaluation of In[91]:= NIntegrate::nlim: etaPrime = \[Eta] is not a valid limit of integration.

During evaluation of In[91]:= FindRoot::nlnum: The function value {-1.+NIntegrate[scatterRate[etaPrime,\[Rho]],{etaPrime,\[Eta],eta0}]} is not a list of numbers with dimensions {1} at {\[Eta]} = {5.17952*10^16}.

During evaluation of In[91]:= ReplaceAll::reps: {FindRoot[\[Tau][\[Eta],\[Rho]]==1,{\[Eta],etaZ[1100]}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

During evaluation of In[91]:= NIntegrate::nlim: etaPrime = \[Eta] is not a valid limit of integration.

During evaluation of In[91]:= General::stop: Further output of NIntegrate::nlim will be suppressed during this calculation.

During evaluation of In[91]:= FindRoot::nlnum: The function value {-1.+NIntegrate[scatterRate[etaPrime,\[Rho]],{etaPrime,\[Eta],eta0}]} is not a list of numbers with dimensions {1} at {\[Eta]} = {5.17952*10^16}.

During evaluation of In[91]:= General::stop: Further output of FindRoot::nlnum will be suppressed during this calculation.

During evaluation of In[91]:= ReplaceAll::reps: {FindRoot[\[Tau][\[Eta],\[Rho]]==1,{\[Eta],etaZ[1100]}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

During evaluation of In[91]:= General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation.

The evaluation of the function looks fine. I can put in any reasonable number and the output is exactly what is expected. The plot, also, looks fine in the given range and gives me the solution I'm searching for if I use a ruler: enter image description here

So how come, when I try to evaluate a 'FindRoot' or 'NMinimize', I get this mass of errors. I can't make sense of them and am especially frustrated because the plot and single execution of the function don't produce the errors. Also, the function runs for hours and gives me back the starting value. What am I doing wrong?

$\endgroup$
3
  • $\begingroup$ Use Block in AngularScale for the local read-only vairable z. AngularScale[ρB_] := Block[{z = zEta[etaStar[ρB]]},sStar[z, ρB]/DLSS[z]]. The errors appear to be caused by evaluation order / using symbols in NIntegrate in the lines τ[η_, ρB_] := NIntegrate[scatterRate[etaPrime, ρB], {etaPrime, η, eta0}] etaStar[ρB_] := η /. Quiet[FindRoot[τ[η, ρB] == 1, {η, etaZ[1100]}]] It may help to change some of your functions like tau so they only accept NumericQ arguments. $\endgroup$
    – flinty
    Jun 14 '20 at 20:15
  • $\begingroup$ @flinty - Yes, I agree the error appears in the integration function, but why does it appear when I perform optimization and not when I plot or calculate individual values? btw - the 'Block' appears to have cleared some of the errors out, but not the big ones. $\endgroup$
    – Quarkly
    Jun 14 '20 at 20:28
  • $\begingroup$ Because of this I think mathematica.stackexchange.com/a/26037/72682 $\endgroup$
    – flinty
    Jun 14 '20 at 20:36
4
$\begingroup$

I've got: {1.03051, {ρ -> 5.*10^-27}} after 15 minutes or so - I lost track of time. I got it by using Block in AngularScale as mentioned in the comments, but most importantly by plastering your code with NumericQ in lots of places where you have arguments:

Clear["Global`*"]
c = 2.99792458*^8;
A = 4.38*^-11;
t0 = 8.57697*10^17;
mB = 2.91*^-27;
ργ = 4.64*^-31;
σT = 6.65*^-29;
eta0 = 1.82*^18;

FreeElectronFractionData = {{3000, 1.0829044`}, {2909.991`, 
    1.0827902`}, {2819.9821`, 1.0825628`}, {2729.9731`, 
    1.0818148`}, {2639.9641`, 1.0792593`}, {2549.9551`, 
    1.0709482`}, {2459.9462`, 1.0503868`}, {2369.9372`, 
    1.0224601`}, {2279.9282`, 1.0056397`}, {2189.9193`, 
    1.0009301`}, {2099.9103`, 1.0001211`}, {2009.9013`, 
    1.0000111`}, {1919.8923`, 0.9999931`}, {1829.8834`, 
    0.99996694`}, {1739.8744`, 0.99984236`}, {1649.8654`, 
    0.99910856`}, {1559.8565`, 0.99387158`}, {1469.8475`, 
    0.95045371`}, {1379.8385`, 0.69166608`}, {1289.8296`, 
    0.2628594`}, {1199.8206`, 0.059562751`}, {1109.8116`, 
    0.0097843423`}, {1019.8026`, 0.0013015781`}, {929.79367`, 
    0.00017836243`}, {839.7847`, 0.000036054186`}, {749.77573`, 
    0.000012683805`}, {659.76676`, 7.0873662`*^-6}, {569.75779`, 
    5.0157033`*^-6}, {479.74881`, 3.98173`*^-6}, {389.73984`, 
    3.377567`*^-6}, {299.73087`, 2.9942744`*^-6}, {209.7219`, 
    2.740874`*^-6}, {119.71293`, 2.5702004`*^-6}, {29.70396`, 
    2.4566939`*^-6}};

etaT[t_?NumericQ] := t*(2*c + A*t)/c
tZ[z_?NumericQ] := Sqrt[t0^2/(z + 1)]
etaZ[z_?NumericQ] := etaT[tZ[z]]
zEta[η_?NumericQ] := 
 t0^2/((-c + Sqrt[c]*Sqrt[c + A*η])/A)^2 - 1
a[z_?NumericQ] := 1/(z + 1)

R[z_?NumericQ, ρB_?
   NumericQ] := (3*ρB)/((a[z]^3*(4*ργ))/a[z]^4)
freeElectronFraction = Interpolation[FreeElectronFractionData];
electronDensity[
  z_?NumericQ, ρB_?
   NumericQ] := (freeElectronFraction[z]*ρB*(1 + z)^3)/mB
scatterRate[η_?NumericQ, ρB_?NumericQ] := 
 electronDensity[zEta[η], ρB]*σT*a[zEta[η]]*c
τ[η_?NumericQ, ρB_?NumericQ] := 
 NIntegrate[scatterRate[etaPrime, ρB], {etaPrime, η, eta0}]
etaStar[ρB_?NumericQ] := η /. 
  Quiet[FindRoot[τ[η, ρB] == 1, {η, etaZ[1100]}]]
sStar[z_?NumericQ, ρB_?NumericQ] := 
  NIntegrate[
   c/Sqrt[3*(1 + R[zEta[etaPrime], ρB])], {etaPrime, 0, etaZ[z]}];
DLSS[z_?NumericQ] := (t0*(A*t0*z + 2*c*(1 + z - Sqrt[1 + z])))/(2 + z)
AngularScale[ρB_?NumericQ] := Block[{z = zEta[etaStar[ρB]]},
  sStar[z, ρB]/DLSS[z]]

AngularScale[5.453*10^-27]


Plot[AngularScale[ρ], {ρ, 5*10^-27, 6*10^-27}]

(* Note - Modified this from the original *)
NMinimize[{Abs[
   AngularScale[ρ] - 0.010411], ρ > 5*10^-27 && ρ < 
    6*10^-27}, {ρ}]

$\endgroup$
6
  • $\begingroup$ That's it. Thank you so much. I can't imagine, yet, what kind of symbols the FindMinimum, NMinimize and FindRoot functions are jamming into my functions, but I guess now at least I know what part of the documentation I need to dive into. $\endgroup$
    – Quarkly
    Jun 14 '20 at 21:14
  • $\begingroup$ Yes I hate it but NumericQ is required so often: support.wolfram.com/12502 When in doubt, if NIntegrate/FindRoot/NMinimize etc. are failing due to symbols, it's probably some argument x to a function that needs to be x_?NumericQ $\endgroup$
    – flinty
    Jun 14 '20 at 21:18
  • $\begingroup$ It took me a few minutes to find this, but you changed the 'Integrate' on sStar to 'NIntegrate'. That seems to have made the biggest difference in the functionality of this thing. Why didn't 'Integrate' work in this situation? $\endgroup$
    – Quarkly
    Jun 14 '20 at 22:04
  • $\begingroup$ Oh yes, Integrate was trying to find a symbolic antiderivative - I forgot to mention that. $\endgroup$
    – flinty
    Jun 14 '20 at 22:20
  • $\begingroup$ I never would have caught that. $\endgroup$
    – Quarkly
    Jun 14 '20 at 22:38

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.