# Solving nested integral equations with NSolve

I want to solve the following integral equations for $$m$$ and $$q$$, with parameters $$J_0$$ and $$T$$: $$q(J_0,T)=\int_{-\infty}^{\infty} \frac{1}{2\pi}e^{-\frac{z^2}{2}}\,\mbox{tanh}^2\left(\frac{z\sqrt{q(J_0,T)}+J_0 m(J_0,T)}{T}\right) dz$$

$$m(J_0,T)=\int_{-\infty}^{\infty} \frac{1}{2\pi}e^{-\frac{z^2}{2}}\,\mbox{tanh}\left(\frac{z\sqrt{q(J_0,T)}+J_0 m(J_0,T)}{T}\right) dz$$

There is a way to solve them using FindRoot, but the computation time is very large, so I need an alternative. What I wanted to use is NSolve which can solve equations for multiple variables at the same time. I tried with the following code:

precision = 100;
Maxit = 500;
accuracy = Round[precision/2];
NSolve[{
q == NIntegrate[
1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[SetPrecision[(Sqrt[q] z + 0.5 m)/0.5, precision + 1]]^2,
{z, -Infinity, Infinity},
WorkingPrecision -> precision,
AccuracyGoal -> accuracy],
m == NIntegrate[
1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[SetPrecision[(Sqrt[q] z + 0.5 m)/0.5, precision + 1]],
{z, -Infinity, Infinity},
WorkingPrecision -> precision, AccuracyGoal -> accuracy] },
{q, m}]


but it seems that NSolve cannot use NIntegrate since it doesn't give any result. Can someone help me?

• This looks very similar to your other question Commented May 23, 2023 at 15:58
• (1) FindRoot can solve for multiple variables at the same time. (2) NSolve is not designed to solve such problems as this. (3) Why not code the Tanh[] factor as Tanh[2 (Sqrt[q] z + m/2)]? (It won't solve the root problem, but it makes the code easier to manage, imho.) Commented May 23, 2023 at 16:01
• I worked on the original problem. My conclusion was that the plots given in that post were quite accurate. However it's based on using FindRoot which in my view is problematic as the routine may not converge to the desired root depending on what basin the seed is in. However, if instead you used the values of that algorithm as seeds to a second iteration of FindRoot then the convergence problem is greatly reduced.
– josh
Commented May 23, 2023 at 19:52

Somehow I thought this would be harder. Did I make a mistake?

NumericQ-protected integrals:

ClearAll[int, int2];
int[m_?NumericQ, q_?NumericQ] :=
NIntegrate[
1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[2 (Sqrt[q] z + m/2)],
{z, -Infinity, Infinity},
WorkingPrecision -> precision, AccuracyGoal -> accuracy];
int2[m_?NumericQ, q_?NumericQ] :=
NIntegrate[
1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[2 (Sqrt[q] z + m/2)]^2,
{z, -Infinity, Infinity},
WorkingPrecision -> precision, AccuracyGoal -> accuracy];


Low-precision trial:

precision = MachinePrecision;
accuracy = Round[precision/2];
FindRoot[{q == int2[m, q], m == int[m, q]}, {{q, 1}, {m, 1/10}},
WorkingPrecision -> precision]

(*  {q -> 0.530368, m -> -1.17503*10^-16}  *)


High-precision trial, with starting points from low-precision trial (wonder if m == 0?):

precision = 100;
accuracy = Round[precision/2];
FindRoot[{q == int2[m, q], m == int[m, q]}, {{q, 53/100}, {m, 1/100}},
WorkingPrecision -> precision]
(*
{q -> 0.5303683920507946325937163727372733258213296547926817305857681565194767283477281152301036248598861080,
m -> -3.839206009833973918466466394691297613302158620963752726865098462765893379586940948647969930045448362*10^-116}
*)


Try with m == 0:

FindRoot[{q == int2[0, q]}, {{q, 53/100}},
WorkingPrecision -> precision]

(*  {q -> 0.5303683920507946325937163727372733258213296547926817305857681565194767283477281152301036248598861080}  *)


Check:

int[m, q] /. % /. m -> 0

(*  0.*10^-101  *)

• Don't see how that code is giving the variables $J$ and $T$ for the computed values of $q$ and $m$.
– josh
Commented May 23, 2023 at 19:47
• @josh The OP substituted special values for the parameters. If the OP wants a different code to work, they should have posted it instead of what they did, imo. It's sort of a mess, and I'm not wading into that tar pit.... Commented May 23, 2023 at 19:53
• Oh I see that now. Sorry. Beautiful tar bit nevertheless. :)
– josh
Commented May 23, 2023 at 20:20
• I tried with specific values but the final goal is to compute $m$ and $q$ for $0<J_0<2$ and $0<T<2$, as @josh saw in the other post. Actually I needed a more efficient (less time consuming) way to tackle the problem to another similar situation where there are more nested integral equations Commented May 24, 2023 at 7:38

Alternatively to MichaelE2's interesting answer you might use

- NMinimize

J[q_?NumericQ,
m_?NumericQ] := # . # &[{q -
NIntegrate[
1/Sqrt[2 Pi] Exp[-z^2/
2] Tanh[ (Sqrt[q] z + 1/2 m)/(1/2) ]^2, {z, -Infinity, -10,
m - NIntegrate[
1/Sqrt[2 Pi] Exp[-z^2/
2] Tanh[ (Sqrt[q] z + 1/2 m)/(1/2)], {z, -Infinity, -10, 10,

NMinimize[ Re@J[q, m], {q, m}] // Quiet
(*{1.19349*10^-17, {q -> 0.530368, m -> 6.34778*10^-9}}*)


or

- FixedPointList

intqm[q_?NumericQ, m_?NumericQ] := {
NIntegrate[1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[ (Sqrt[q] z + 1/2m)/(1/2)]^2, {z, -Infinity, -10, 10, Infinity}, Method -> "GlobalAdaptive"],
NIntegrate[1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[ (Sqrt[q] z + 1/2m)/(1/2)],{z, -Infinity, -10, 10, Infinity}, Method -> "GlobalAdaptive"]}

FixedPointList[Apply[intqm, #] &, {1, 1}, 15]
(*{..., {0.530368, 0.0000164061}, {0.530368, 7.70484*10^-6}}*)


Both results agree well with MichaelE2's answer!