Need help solving iterativelly a set of transcendental equations with NIntegrate

Question

I have two sets of equations:

$q = \frac{1}{\sqrt{2 \pi}}\int dz~ e^{-z^2/2}F(q,\phi;z,\beta,\Delta, J) $

$\phi = \frac{1}{\sqrt{2 \pi}}\int dz~ e^{-z^2/2}G(q,\phi;z,\beta, \Delta, J) $

As you may see they are self-consistent equations, the main problem is that the expressions for $F$ and $G$ quite involved, so the integration over $z$ can only be achieved by a numerical integration.

My goal is to these self-consistent equations for several values of the parameters $\Delta$ and $\beta$, then plot $x(\beta)$ for different values of $\Delta$. The problems I'm proceeding in the following way:

The functions $F$ and $G$ are actually fractions with a common denominator I define as

Den[\[CapitalDelta]_, J_, \[Beta]_, q_, \[Phi]_, z_] = 
 (1/\[CapitalDelta])*Integrate[Exp[(-2^(-1))*(q - \[Phi])*\[Beta]^2*J^2*t^2]*
    Cosh[\[Beta]*J*Sqrt[q]*z*t], {t, 1 - \[CapitalDelta]/2, 1 + \[CapitalDelta]/2}]

Then we comput the full functions $F$ and $G$, defined as fractionq and fraction$\phi$ in the latex code

fractionq[\[CapitalDelta]_, J_, \[Beta]_, q_, \[Phi]_] = 1/\[CapitalDelta]*Integrate[Exp[-1/2 (q-\[Phi]) \[Beta]^2 J^2 t^2] (t^2 Cosh[\[Beta] J Sqrt[q] z t]- z t /(\[Beta] J Sqrt[q])Sinh[\[Beta] J Sqrt[q]z t]),{t, 1-\[CapitalDelta]/2,1+\[CapitalDelta]/2}]/Den[\[CapitalDelta], J, \[Beta], q, \[Phi], z]
fraction\[Phi][\[CapitalDelta]_, J_, \[Beta]_, q_, \[Phi]_] = 1/\[CapitalDelta]*Integrate[Exp[-1/2 (q-\[Phi]) \[Beta]^2 J^2 t^2] t^2 Cosh[\[Beta] J Sqrt[q] z t],{t, 1-\[CapitalDelta]/2,1+\[CapitalDelta]/2}]/Den[\[CapitalDelta], J, \[Beta], q, \[Phi], z]

We then prepare the HoldForm of the numerical integration over $z$

RHSq[\[CapitalDelta]_, J_, \[Beta]_, q_, \[Phi]_]:= NIntegrate[1/Sqrt[2 Pi] Exp[-1/2 z^2] fractionq[\[CapitalDelta], J, \[Beta], q, \[Phi]], {z, -Infinity, Infinity},
 PrecisionGoal->4, WorkingPrecision->10, AccuracyGoal->10, MinRecursion->2]
RHS\[Phi][\[CapitalDelta]_, J_, \[Beta]_, q_, \[Phi]_]:= NIntegrate[1/Sqrt[2 Pi] Exp[-1/2 z^2] fraction\[Phi][\[CapitalDelta], J, \[Beta], q, \[Phi]], {z, -Infinity, Infinity}, 
PrecisionGoal->4, WorkingPrecision->10, AccuracyGoal->10, MinRecursion->2]

And finally we defined the ordered pair {RHSq, RHS$\phi$} to which we are going to apply a function (either FixedPoint, Nest, NestWhile) in order to solve iteratively for several values of $\beta$ and $\Delta$.

newpairs[\[CapitalDelta]_, J_, \[Beta]_][{q_, \[Phi]_}]:={RHSq[\[CapitalDelta], J, \[Beta], q, \[Phi]], RHS\[Phi][\[CapitalDelta], J, \[Beta], q, \[Phi]]}

Table[FixedPoint[newpairs[0.2, 1, \[Beta]],{1.5,1}],{\[Beta],\[Beta]list}]

where $\beta$list is just a list of the different values of $\beta$ I want to generate.

TemperatureList = Table[i, {i, 0.0002, 2.0002, 0.02}]
\[Beta]list = 1/TemperatureList

The following are the error outputs displayed while running this routine:

NIntegrate::inumri "The integrand ... has evaluated to Overflow, Indeterminate, or Infinity for all
sampling points in the region with boundaries
{{-[Infinity],-3.000000000}}"

General::stop "Further output of NIntegrate::inumri will be suppressed during this
calculation"

NIntegrate::write "tag Times in -z is Protected"

several times

Answer 1

Your code can be made to work by one major change and a couple of smaller changes.

• FindRoot is generally more efficient than FixedPoint for finding roots of equations. Unless FixedPoint is necessary for the particular physical situation you're trying to model (renormalization flows?), I would recommend you use FindRoot instead.

• Many of the errors in your code are caused by calls to NIntegrate with symbolic integrands, rather than integrands that evaluate to a particular number when provided with an input value in the integration range. This can be avoided by defining the input slots with PatternTests:

  fractionq[\[CapitalDelta]_?NumericQ, J_?NumericQ, \[Beta]_?NumericQ, q_?NumericQ, \[Phi]_?NumericQ, z_?NumericQ] 
     := (* long expression including NIntegrate *)
  

  This syntax means that the function (and by extension, NIntegrate) will only be called when all of the arguments provided are numerical expressions; otherwise, it will be returned unevaluated.

• To do this input sanitization properly, the argument z should be an argument to the integrand functions fractionq and fraction\[Phi].

With these modifications, Mathematica returns a result for q and \[Phi] in 20–140 seconds (on my five-year-old Macbook Pro) for values of $\beta$ between 0.5 and 2, with lower $\beta$ taking longer. This would make it feasible to run this code for the 100 data points that you want; your result would be ready in a few hours.

Further speed-ups could be obtained by using ParallelTable (allowing your computer to calculate several data points at once); using an "annealing" process where the results for one value of \[Beta] are used as initial guesses for the next value of \[Beta] (reducing the number of steps FindRoot has to make before getting convergence); and making compiled versions of the integrand functions.

---

Modified code:

fractionq[\[CapitalDelta]_?NumericQ, J_?NumericQ, \[Beta]_?NumericQ, 
  q_?NumericQ, \[Phi]_?NumericQ, z_?NumericQ] := 
 1/\[CapitalDelta]*
  NIntegrate[
    Exp[-1/
        2 (q - \[Phi]) \[Beta]^2 J^2 t^2] (t^2 Cosh[\[Beta] J Sqrt[
           q] z t] - 
       z t/(\[Beta] J Sqrt[q]) Sinh[\[Beta] J Sqrt[q] z t]), {t, 
     1 - \[CapitalDelta]/2, 1 + \[CapitalDelta]/2}]/
   Den[\[CapitalDelta], J, \[Beta], q, \[Phi], z]

fraction\[Phi][\[CapitalDelta]_?NumericQ, 
  J_?NumericQ, \[Beta]_?NumericQ, q_?NumericQ, \[Phi]_?NumericQ, 
  z_?NumericQ] := 
 1/\[CapitalDelta]*
  NIntegrate[
    Exp[-1/2 (q - \[Phi]) \[Beta]^2 J^2 t^2] t^2 Cosh[\[Beta] J Sqrt[
        q] z t], {t, 1 - \[CapitalDelta]/2, 1 + \[CapitalDelta]/2}]/
   Den[\[CapitalDelta], J, \[Beta], q, \[Phi], z]

RHSq[\[CapitalDelta]_?NumericQ, J_?NumericQ, \[Beta]_?NumericQ, 
  q_?NumericQ, \[Phi]_?NumericQ] := 
 NIntegrate[
  1/Sqrt[2 Pi] Exp[-1/2 z^2] fractionq[\[CapitalDelta], J, \[Beta], 
    q, \[Phi], z], {z, -Infinity, Infinity}, PrecisionGoal -> 4, 
  WorkingPrecision -> 10, AccuracyGoal -> 10, MinRecursion -> 2]

RHS\[Phi][\[CapitalDelta]_?NumericQ, J_?NumericQ, \[Beta]_?NumericQ,     
  q_?NumericQ, \[Phi]_?NumericQ] := 
 NIntegrate[
  1/Sqrt[2 Pi] Exp[-1/2 z^2] fraction\[Phi][\[CapitalDelta], 
    J, \[Beta], q, \[Phi], z], {z, -Infinity, Infinity}, 
  PrecisionGoal -> 4, WorkingPrecision -> 10, AccuracyGoal -> 10, 
  MinRecursion -> 2]

newpairs[\[CapitalDelta]_, 
   J_, \[Beta]_][{q_?NumericQ, \[Phi]_?
    NumericQ}] := {RHSq[\[CapitalDelta], J, \[Beta], q, \[Phi]], 
  RHS\[Phi][\[CapitalDelta], J, \[Beta], q, \[Phi]]}

\[Beta] = 1.;
FindRoot[newpairs[0.2, 1, \[Beta]][{q, \[Phi]}] == {q, \[Phi]}, {{q, 1}, {\[Phi], 1.5}}]

(* {q -> 0.0101124, \[Phi] -> 1.01008} *)