1
$\begingroup$

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?

$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – josh
    May 4, 2023 at 17:26
  • $\begingroup$ 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! $\endgroup$ May 4, 2023 at 20:06
  • $\begingroup$ 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. $\endgroup$
    – josh
    May 4, 2023 at 20:57
  • $\begingroup$ There is a typo with $e^{z^2/2}$ in your Latex expressions, it should be $e^{-z^2/2}$. $\endgroup$ May 5, 2023 at 7:26
  • $\begingroup$ Thanks for that, I'll correct it $\endgroup$ May 5, 2023 at 7:42

1 Answer 1

3
$\begingroup$

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]

![enter image description here

Hope it helps!

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.