# 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 May 23 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.) May 23 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
May 23 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
May 23 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.... May 23 at 19:53
• Oh I see that now. Sorry. Beautiful tar bit nevertheless. :)
– josh
May 23 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 May 24 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!