Why does NMinimize start spouting error messages?

Question

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:

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?

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. — flinty, Jun 14 '20 at 20:15
@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. — Quark Soup, Jun 14 '20 at 20:28
Because of this I think https://mathematica.stackexchange.com/a/26037/72682 — flinty, Jun 14 '20 at 20:36

score 4 · Accepted Answer · edited Jun 14 '20 at 21:19

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}, {ρ}]

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. — Quark Soup, Jun 14 '20 at 21:14
Yes I hate it but NumericQ is required so often: https://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 — flinty, Jun 14 '20 at 21:18
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? — Quark Soup, Jun 14 '20 at 22:04
Oh yes, Integrate was trying to find a symbolic antiderivative - I forgot to mention that. — flinty, Jun 14 '20 at 22:20
Do you know what's odd? FindMinimum and FindRoot appear to choke on the same formula that NMinimize is using. Any idea why FindRoot doesn't work here? — Quark Soup, Jun 14 '20 at 22:41

Why does NMinimize start spouting error messages?

1 Answers1