2
$\begingroup$

I have a function SQ[b,zQ,zh] where I want to find at which zQ it is a minimum given b and zh. I set b=10^-x and zh=1.355 10^-x where I choose x. For x=3 I encounter a FindMinimum::lstol issue. I have tried changing MaxRecursion and WorkingPrecision but it is not resolved. Can anyone help me with this?

I also want to extend the calculation for larger x=4,5,...,10

d = 3;
ag = 10;
pg = 10;
wp = 50;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
torootsig[b_?NumericQ, sig_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := Module[{br, sigr, zQr, zhr}, {br, sigr, zQr, zhr} = Rationalize[{b, sig, zQ, zh}, 0]; br - NIntegrate[z^d/Sqrt[f[z, zhr] (zQr^(2 d) (1 + (sigr^2/f[zQr, zhr])) - z^(2 d))], {z, 0, zQr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 100]]
sig[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := sig /. FindRoot[torootsig[b, sig, zQ, zh] == 0, {sig, -50 1.3 10^-3, 0}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 100]
intSQ1[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := Module[{br, zQr, zhr}, {br, zQr, zhr} = Rationalize[{b, zQ, zh}, 0]; (-1/(d - 1)) (1/(zQr^(2 d) (1 + sig[br, zQr, zhr]^2/f[zQr, zhr]))) NIntegrate[z^d Sqrt[f[z, zhr]/(1 - (1/(zQr^(2 d) (1 + sig[br, zQr, zhr]^2/f[zQr, zhr]))) z^(2 d))], {z, 0, zQr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 100]]
intSQ2[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := Module[{br, zQr, zhr}, {br, zQr, zhr} = Rationalize[{b, zQ, zh}, 0]; (-1/(2 zhr^(d + 1))) ((d + 1)/(d - 1)) NIntegrate[z Sqrt[(1 - (1/(zQr^(2 d) (1 + sig[br, zQr, zhr]^2/f[zQr, zhr]))) z^(2 d))/f[z, zhr]], {z, 0, zQr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 100]]
intSQ3[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := Module[{br, zQr, zhr}, {br, zQr, zhr} = Rationalize[{b, zQ, zh}, 0]; (1/zhr)^(d + 1) NIntegrate[z/Sqrt[f[z, zhr] (1 - (1/(zQr^(2 d) (1 + sig[br, zQr, zhr]^2/f[zQr, zhr]))) z^(2 d))], {z, 0, zQr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 300, Method -> "LocalAdaptive"]]
SQ[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := Module[{br, zQr, zhr}, {br, zQr, zhr} = Rationalize[{b, zQ, zh}, 0]; (-Sqrt[f[zQr, zhr] (1 - (1/(zQr^(2 d) (1 + sig[br, zQr, zhr]^2/f[zQr, zhr]))) zQr^(2 d))]/((d - 1) zQr^(d - 1)) + intSQ1[br, zQr, zhr] + intSQ2[br, zQr, zhr] + intSQ3[br, zQr, zhr] + 1/zQr^(d - 1))/4 ]

x = 3;
FindMinimum[{SQ[10^-x, zQ, 1.355 10^-x]/10^(2 x)}, {zQ, 0.9955 1.355 10^-x, 0.995 1.355 10^-x, 0.996 1.355 10^-x}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 500] // AbsoluteTiming
FindMinimum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than 50.` digits of working precision to meet these tolerances.
{598.3201164, {0.19195144250523353513414807356587315266329459419565, {zQ -> 0.0013490146179865683928983415239627338448699548869898}}}

There is definitely a minimum as is shown below,

x = 3;
Plot[SQ[10^-x, n 1.355 10^-x, 1.355 10^-x]/10^(2 x), {n, 0.995, 0.996}, PlotStyle -> {Blue, Thickness[0.005]}, PlotRange -> Full, ImageSize -> Large] // AbsoluteTiming

Image

Improve this question
$\endgroup$
2
  • $\begingroup$ Your functions are numerical therefore the optimal solution can be computed with NMinimize. $\endgroup$
    – Alex Trounev
    Nov 16, 2021 at 6:00
  • $\begingroup$ @AlexTrounev I have tried NMinimize with the same range given above (excluding initial point) and it just keeps running for so long and not ending so I just aborted it around 30 minutes later. That's so exaggerated for only one point. $\endgroup$
    – mathemania
    Nov 16, 2021 at 8:02

2 Answers 2

Reset to default
4
+50
$\begingroup$

We can simplifier this code to avoid recursive usage of some functions during evaluation as follows

d = 3;
pg = 10;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
torootsig[b_?NumericQ, sig_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := 
  Module[{br = b, sigr = sig, zQr = zQ, zhr = zh}, 
   i1 = br - 
     NIntegrate[
      z^d/Sqrt[
        f[z, zhr] (zQr^(2 d) (1 + (sigr^2/f[zQr, zhr])) - 
           z^(2 d))], {z, 0, zQr}, PrecisionGoal -> pg]; i1];
sig[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := 
  sig /. FindRoot[
    torootsig[b, sig, zQ, zh] == 0, {sig, -50 1.3 10^-3, 0}];
intSQ1[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := 
  Module[{br = b, zQr = zQ, zhr = zh}, 
   i2 = (-1/(d - 
         1)) (1/(zQr^(2 d) (1 + 
           sig[br, zQr, zhr]^2/f[zQr, zhr]))) NIntegrate[
      z^d Sqrt[
        f[z, zhr]/(1 - (1/(zQr^(2 d) (1 + 
                  sig[br, zQr, zhr]^2/f[zQr, zhr]))) z^(2 d))], {z, 0,
        zQr}, PrecisionGoal -> pg]; i2];
intSQ2[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := 
  Module[{br = b, zQr = zQ, zhr = zh}, 
   i3 = (-1/(2 zhr^(d + 1))) ((d + 1)/(d - 1)) NIntegrate[
      z Sqrt[(1 - (1/(zQr^(2 d) (1 + 
                  sig[br, zQr, zhr]^2/f[zQr, zhr]))) z^(2 d))/
         f[z, zhr]], {z, 0, zQr}, PrecisionGoal -> pg]];
intSQ3[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := 
 Module[{br = b, zQr = zQ, zhr = zh}, 
  i4 = (1/zhr)^(d + 1) NIntegrate[
     z/Sqrt[f[z, 
         zhr] (1 - (1/(zQr^(2 d) (1 + 
                 sig[br, zQr, zhr]^2/f[zQr, zhr]))) z^(2 d))], {z, 0, 
      zQr}, PrecisionGoal -> pg, Method -> "LocalAdaptive"]; i4]
SQ[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := 
  Module[{br = b, zQr = zQ, zhr = zh}, 
   i5 = (-Sqrt[
          f[zQr, zhr] (1 - (1/(zQr^(2 d) (1 + 
                    sig[br, zQr, zhr]^2/
                    f[zQr, zhr]))) zQr^(2 d))]/((d - 1) zQr^(d - 1)) +
        intSQ1[br, zQr, zhr] + intSQ2[br, zQr, zhr] + 
       intSQ3[br, zQr, zhr] + 1/zQr^(d - 1))/4; i5];

With this code we can compute list and use Interpolation to plot SQ

x = 3;
lst = Table[{n, 
    10^(-2 x) SQ[10^-x, n 1.355 10^-x, 1.355 10^-x] // Quiet // 
     Chop}, {n, 0.995, 0.996, .0001}];

f0 = Interpolation[lst];

Plot[f0[n], {n, 0.995, 0.996}]

Figure1 Finally we compute min value

FindMinimum[f0[n], {n, .9956}]

Out[]= {0.191951, {n -> 0.995583}}
Improve this answer
$\endgroup$
1
  • $\begingroup$ You used Table to estimate the values near the minimum and used Interpolation to get an estimated curve. From there, you used the interpolated curve to find the minimum. Nice approach, seems accurate. $\endgroup$
    – mathemania
    Nov 17, 2021 at 10:06
4
$\begingroup$

Without defining any precision and accuracy goals and also without pre-setting any working precision we have:

d=3;
f[z_,zh_]:=1-(z/zh)^(d+1);
torootsig[b_?NumericQ,sig_?NumericQ,zQ_?NumericQ,zh_?NumericQ]:=b-NIntegrate[z^d/Sqrt[f[z,zh] (zQ^(2 d) (1+(sig^2/f[zQ,zh]))-z^(2 d))],{z,0,zQ}]
sig[b_?NumericQ,zQ_?NumericQ,zh_?NumericQ]:=sig/.FindRoot[torootsig[b,sig,zQ,zh]==0,{sig,-50 1.3 10^-3,0}]
intSQ1[b_?NumericQ,zQ_?NumericQ,zh_?NumericQ]:=(-1/(d-1)) (1/(zQ^(2 d) (1+sig[b,zQ,zh]^2/f[zQ,zh]))) NIntegrate[z^d Sqrt[f[z,zh]/(1-(1/(zQ^(2 d) (1+sig[b,zQ,zh]^2/f[zQ,zh]))) z^(2 d))],{z,0,zQ},MaxRecursion->100]
intSQ2[b_?NumericQ,zQ_?NumericQ,zh_?NumericQ]:=(-1/(2 zh^(d+1))) ((d+1)/(d-1)) NIntegrate[z Sqrt[(1-(1/(zQ^(2 d) (1+sig[b,zQ,zh]^2/f[zQ,zh]))) z^(2 d))/f[z,zh]],{z,0,zQ},MaxRecursion->100]
intSQ3[b_?NumericQ,zQ_?NumericQ,zh_?NumericQ]:=(1/zh)^(d+1) NIntegrate[z/Sqrt[f[z,zh] (1-(1/(zQ^(2 d) (1+sig[b,zQ,zh]^2/f[zQ,zh]))) z^(2 d))],{z,0,zQ},MaxRecursion->300]
SQ[b_?NumericQ,zQ_?NumericQ,zh_?NumericQ]:=(-Sqrt[f[zQ,zh] (1-(1/(zQ^(2 d) (1+sig[b,zQ,zh]^2/f[zQ,zh]))) zQ^(2 d))]/((d-1) zQ^(d-1))+intSQ1[b,zQ,zh]+intSQ2[b,zQ,zh]+intSQ3[b,zQ,zh]+1/zQ^(d-1))/4

Then we can estimate the minimum for $x=3,4,\ldots,10$ (with a lot of warnings):

 Table[{x,FindMinimum[{SQ[10^-x, zQ, 1.355 10^-x]/10^(2 x)}, {zQ, 0.9955 1.355 10^-x, 0.995 1.355 10^-x, 0.996 1.355 10^-x}]}, {x,Range[3, 10]}]

{{3, {0.191951, {zQ -> 0.00134901}}},{4, {0.191951, {zQ -> 0.000134897}}}, {5, {0.191952, {zQ -> 0.000013489}}}, {6, {0.191952, {zQ -> 1.3489*10^-6}}}, {7, {0.191952, {zQ -> 1.3489*10^-7}}}, {8, {0.191952, {zQ ->1.3489*10^-8}}}, {9, {0.191952, {zQ -> 1.3489*10^-9}}}, {10, {0.191952, {zQ -> 1.3489*10^-10}}}}

I have checked graphically for x=3 and x=4 and the minimum values seem to make sense.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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