NIntegrate inside FindRoot has evaluated to non-numerical values

Question

I have a function Fdiff[m] that I want to put in FindRoot. However, it has given me an error saying that the integrand has evaluated to non-numerical values. I'm aware of other posts having this issue as well, but reading those did not lead to a solution in my case.

I have used ?NumericQ as well as trying to integrate first to get an expression just in terms of m before putting in FindRoot, but I think this can't be done since my integrand is a bit complicated so maybe there's no analytic result, so using NIntegrate is must, but for this to work m has to have a value. So I checked that Fdiff[m] for some specific value of m indeed works.

I believe the syntax of my code is correct but I'm not sure why this problem persist. Any help?

Clear["Global`*"]
a = Rationalize[4.046];
b = Rationalize[0.01613];
c = Rationalize[0.227];
A[z_] := -a Log[b z^2 + 1]
int1[zh_] = ((Integrate[y^3 Exp[-3 A[y]], y] /. {y -> zh}) - (Integrate[y^3 Exp[-3 A[y]], y] /. {y -> 0})) // Simplify;
int2[zh_] = ((Integrate[y^3 Exp[-3 A[y] + c y^2], y] /. {y -> zh}) - (Integrate[y^3 Exp[-3 A[y] + c y^2], y] /. {y -> 0})) // Simplify;

s[zh_] := Exp[3 A[zh]]/(4 zh^3);
T[m_, zh_] := Rationalize[((zh^3 Exp[-3 A[zh]])/(4 Pi int1[zh])) (1 - ((2 c m^2)/(1 - Exp[c zh^2])^2) (Exp[c zh^2] int1[zh] - int2[zh]))] // Simplify;
Tp[m_, zh_] := D[T[m, zh], zh] // Simplify;

Fdiff[m_?NumericQ] := NIntegrate[s[y] Tp[m, y], {y, 3, 100}] - NIntegrate[s[y] Tp[m, y], {y, 7, 100}]

FindRoot[Chop[Fdiff[m]] == 0, {m, 0.048}, AccuracyGoal -> 10, WorkingPrecision -> 10]

NIntegrate::inumr: The integrand (-((y^6 (1+<<1>>)^(6069/<<3>>) <<1>> (1-(227 <<1>> (Power[<<2>>] <<1>>+<<1>>))/(500 (-1+<<1>>)^2)))/(4000000000000000000000000000000000000000000000000000000000000 <<2>> <<1>>^2))+<<1>><<1>><<1>>+<<1>>/<<1>>)/(4 y^3 (1+(1613 y^2)/100000)^(6069/500)) has evaluated to non-numerical values for all sampling points in the region with boundaries {{4,28.}}.

Answer 1

$Version

(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)

Clear["Global`*"]

a = Rationalize[4.046];
b = Rationalize[0.01613];
c = Rationalize[0.227];
A[z_] = -a  Log[b  z^2 + 1];
int1[zh_] = Integrate[y^3  Exp[-3  A[y]], {y, 0, zh}] // Simplify;
int2[zh_] = Integrate[y^3  Exp[-3  A[y] + c  y^2], {y, 0, zh}] // Simplify;

s[zh_] = Exp[3  A[zh]]/(4  zh^3);
T[mv_, zhv_] =
  ((zh^3  Exp[-3  A[zh]])/(4  Pi  int1[
         zh]))  (1 - ((2  c  m^2)/(1 - Exp[c  zh^2])^2)  (Exp[c  zh^2]  int1[
           zh] - int2[zh])) // Simplify;

In the definition of Tp you cannot use D[T[m, zh], zh] when zh is being given a numeric value. Either use Set to evaluate the derivative before zh has a numeric value, or use the form Derivative[0, 1][T][m, zh].

Tp[m_, zh_] = D[T[m, zh], zh] // Simplify;

Fdiff[mv_?NumericQ] := Module[
  {m = Rationalize[mv, 0]},
  NIntegrate[s[y]  Tp[m, y], {y, 3, 7},
   WorkingPrecision -> 15]]

Fdiff /@ {-0.048, 0, 0.048}

(* {-8.33896317973750*10^-6, -8.02294953359420*10^-6, -8.33896317973750*10^-6} *)

Fdiff does not have a root, it is negative

Plot[Fdiff[m], {m, -0.1, 0.1}]

FindMaximum[Fdiff[m], {m, 0}, WorkingPrecision -> 15]

(* {-8.02294953359420*10^-6, {m -> 0}} *)