# Numerical solution of simple nested integral equations

I want to solve the following two equations:

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

and

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

and plot T as a function of h (or h as a function of T). Basically the first equation gives the function $$q(h,T)$$ for each value of $$h$$ and $$T$$ and plugging this result into the second equation we have a relation between $$h$$ and $$T$$.

I tried numerically with the following code

fun[q_?NumericQ, T_?NumericQ, h_?NumericQ] :=  NIntegrate[ 1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[(Sqrt[q] z + h)/T]^2, {z, -Infinity, Infinity}];

qfun[T_?NumericQ, h_?NumericQ] := (FindRoot[q - fun[q, T, h], {q, 1}] // Chop )[[1]][[2]]



which gives the desired results for the function $$q(h,T)$$. Indeed I tried to plot this function getting the expected result by

Plot[{Evaluate@Table[qfun[T, h], {h, {0, 0.1, 0.25, 1.5}}]}, {T, 0.00000001, 2}]


Similarly I tried to numerically integrate and solve the second equation

fun2[h_?NumericQ, T_?NumericQ] := NIntegrate[(1/Sqrt[2 Pi]Exp[-z^2/2] Sech[SetPrecision[(Sqrt[qfun[T, h]] z + h)/T, 31]]^4) // Rationalize, {z, -Infinity, Infinity}, AccuracyGoal -> 20, PrecisionGoal -> 20, WorkingPrecision -> 30];

hfun[T_?NumericQ] := (FindRoot[T^2 - fun2[h, T] , {h, 0.4}] // Chop)[[1]][[2]]


But then, plotting $T$ as a function of $h$

Plot[Tfun[h], {h, 0.01, 1.5}]


I didn't get the expected result. In particular the plot is:

which is quite similar to what I expect, but it seems that for some values of $$h$$ the function FindRoot gives the right value for $$T$$ but with negative sign.

I tried with Chop and with some changes of numerical precision but without success, can anyone help me?

• When I use tfun[h_?NumericQ] := (FindRoot[T^2 - fun2[h, T], {T, 0.4}] Chop)[[1]][[2]]; hPoints = Table[{h, tfun[h]}, {h, 0.02, 1.5, 0.1}]; then ListPlot[hPoints] appears continuous like the upper curve in your plot. Is this the form of T(h) you're using?
– josh
Commented May 4, 2023 at 17:26
• Yes, the expected behaviour of T(h) is the one resulting from the code you suggested, thank you!. By the way I was wondering if there is a way to systematically get the "continuous" result, because then I want to compute more complicated things. Thank you anyway! Commented May 4, 2023 at 20:06
• If you're willing to settle for an approximate fit of the data, then can use theF=FindFormula[hpoints,x,n] to find n fits for the data. One which fit well was $1.04268\, -0.719036 \sqrt{x}$. Type FindFormula in the notebook, hover over it, and select the "I" for infor about FindFormula and other fit functions.
– josh
Commented May 4, 2023 at 20:57
• There is a typo with $e^{z^2/2}$ in your Latex expressions, it should be $e^{-z^2/2}$. Commented May 5, 2023 at 7:26
• Thanks for that, I'll correct it Commented May 5, 2023 at 7:42

You have two equations and three unknowns {q,T,h}, though a one dimensional solution {q,T[q],h[q]} is to be expected.

eqn[q_?NumericQ, T_?NumericQ, h_?NumericQ] :={
q -NIntegrate[1/Sqrt[2 Pi] Exp[-z^2/2] Tanh[(Sqrt[q] z + h)/T]^2, {z,-Infinity, -7, 0, 7,Infinity}, PrecisionGoal -> 4, AccuracyGoal -> 6]
,
T^2 - NIntegrate[1/Sqrt[2 Pi] Exp[-z^2/2] Sech[(Sqrt[q] z + h)/T]^4, {z,-Infinity, -7, 0, 7, Infinity}, PrecisionGoal -> 4, AccuracyGoal -> 6]
}


We solve the equations for given parameter 0<qq<1

sol[qq_] := {qq, T, h, # . # &[eqn[qq, T, h]]} /. FindRoot[eqn[qq, T, h] == 0, {T, 1, 0, 2}, {h, 1, 0, 2}]  (*q,T,h,res*)
qth = Table[sol[qq], {qq, Range[0, 1 , #] &[.05]}];


If last entry of qth is zero, FindRoot found a solution!

Graphics3D[{punkte = Select[qth, #[[-1]] < 10^-7 &][[All, {2, 3, 1}]];
Point[punkte], Line[punkte]}, AxesLabel -> {T, h, q}, Axes -> True]


Hope it helps!

• Apart from the exponent of the Sech function, which is 4, not 2, it gives the expected result! Thank you! Commented May 5, 2023 at 15:52
• @saveriopalazzi Thanks for the hint, I corrected my answer Commented May 5, 2023 at 20:28