# Need help solving iterativelly a set of transcendental equations with NIntegrate

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

• Welcome to Mathematica! Please edit your question to include the code you've come up with so far; it will be much easier for users to give you useful advice if they know the details of your problem, particularly the functions $F$ and $G$. Commented Jun 9, 2022 at 18:25
• Thanks for the advice! Sorry this is my first post, do you think is better in this way or can I make it more readable? Commented Jun 10, 2022 at 9:24
• That's much more answerable now, but it throws an error on my system because it doesn't include a definition of βlist. Can you edit your question to include it? Commented Jun 10, 2022 at 12:47
• Yeah! sorry, it's done Commented Jun 10, 2022 at 14:36
• Is there a particular reason you need to use FixedPoint rather than FindRoot? On my five-year-old MacBook, FindRoot returns a root of the equation newpairs[0.2, 1, 1][{q, \[Phi]}] == {q, \[Phi]} (i.e., $\beta = 1$) in a little over one minute, which seems like a reasonable time frame if you want to do 100 $\beta$ values. FixedPoint takes longer (it's been running for five minutes so far without an output), but you may have some particular reason that you need to use it. Commented Jun 10, 2022 at 17:09

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} *)

• Thanks a lot for the help! it seems to be working way better, just that for some values of $\beta$, for example 5000 or even 50, it outputs the following error Wolfram NIntegrate::inumri: The integrand 0.398942 E^(-0.5 z^2) fractionq[0.2,1.,5000.,0.,0.,z] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{-\[Infinity],-3.000000000}}.  I imagine I can correct this error by playing with the NIntegrate parameters right?, which parameters are better to tune in this case? Commented Jun 13, 2022 at 13:31
• @GreivinAlfaro: That I'm not sure about. I would take a look at what the integrand is actually doing in those cases and perhaps look at truncating the range of integration based on that. Also, are you still using FixedPoint for some reason? I only got that error when I tried using FixedPoint, as in your original code. Using FindRoot, as in my modified code, doesn't throw that error at all. Commented Jun 13, 2022 at 14:21
• My guess is that the problem comes from the fact that your integrands are of the form $e^{-\beta^2 k} \cosh (\beta z l)$, where $k$ and $l$ are combinations of other constants. These could easily evaluate to something of the form $0 \times \infty$ internally when Mathematica plugs in real numbers. Commented Jun 13, 2022 at 14:27