1
$\begingroup$

I have a function SQ[b,zQ,zh] where I want to find at which zQ it is a minimum. The parameters are b=0.1, zh=(0.1,0.9) increment by 0.1 so zh has 9 values.

d = 3;
ag = 10;
pg = 10;
wp = 20;
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 -> 20]]
sig[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := sig /. FindRoot[torootsig[b, sig, zQ, zh], {sig, -2 zh, 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 -> 20]]
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 -> 20]]
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 -> 20]]
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 ]

FindMinimum[SQ[0.1, zQ, 0.1], {zQ, 0.0993, 0.0999}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
{41.70848808, {zQ -> 0.09987512206}}

0.09987512206/0.1
0.9987512206

FindMinimum[SQ[0.1, zQ, 0.2], {zQ, 0.195, 0.199}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
{7.388269649, {zQ -> 0.197303209}}

0.197303209/0.2
0.986516045

FindMinimum[SQ[0.1, zQ, 0.3], {zQ, 0.29, 0.299}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
{2.875769877, {zQ -> 0.291574518}}

0.291574518/0.3
0.97191506

FindMinimum[SQ[0.1, zQ, 0.4], {zQ, 0.35, 0.4}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
{1.516662443, {zQ -> 0.3843298312}}

0.3843298312/0.4
0.960824578

FindMinimum[SQ[0.1, zQ, 0.5], {zQ, 0.45, 0.5}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
{0.9366961415, {zQ -> 0.4765573324}}

0.4765573324/0.5
0.9531146648

FindMinimum[SQ[0.1, zQ, 0.9], {zQ, 0.83, 0.86}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
{0.2745876051, {zQ -> 0.8454219108}}

0.8454219108/0.9
0.9393576787

Table[FindMinimum[SQ[0.1,zQ,n/10],{zQ,(93/100) (n/10),(999/1000) (n/10)},AccuracyGoal->ag,PrecisionGoal->pg,WorkingPrecision->wp,MaxIterations->100],{n,1,9,1}]//AbsoluteTiming

I have found at which zQ is the minimum for the two extreme values of zh. For zh=0.1 the minimum is at around zQ ~ 0.998 zh, for zh=0.9 the minimum is at around zQ ~ 0.939 zh.

The problem is as zh goes to lower values FindMinimum will fail for some zQ and below. For example, for zh=0.1, the range I put in FindMinimum is {zQ, 0.0993, 0.0999} but if I change the zQmin=0.0993 to zQmin=0.098, SQ will fail there.

My question is, is there a way to write a code to increment zQmin from some value to some higher value (in this case zQmin = {(93/100) (n/10), (994/1000) (n/10)}) so that when FindMinimum fails it will increment until FindMinimum does not fail. I'm not sure if I'm framing this problem correctly, I was also thinking of maybe Conditionals like If, For, While. My end goal is to execute the Table written above.

My related post will be useful, FindRoot with a big range.

$\endgroup$

1 Answer 1

1
$\begingroup$

In the line Table[FindMinimum[SQ[0.1,zQ,n/10],{zQ,(93/100) (n/10),(999/1000) (n/10)}],{n,1,9,1}] only cases n=1,2 are problematic due to FindRoot since there are no roots for some parameters. For n>=3 we have (all options with ag, pg, wp have been omitted)

d = 3;
ag = 10;
pg = 10;
wp = 20;
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}]];
sig[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := 
  x /. FindRoot[torootsig[b, x, zQ, zh], {x, -2 zh, 0}, 
    MaxIterations -> 1000];
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}]];
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}]];
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}]];
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]
Table[
  FindMinimum[
   SQ[0.1, zQ, n/10], {zQ, (93/100) (n/10), (999/1000) (n/10)}], {n, 
   3, 9, 1}] // AbsoluteTiming

Out[]= {19.1095, {{2.87577, {zQ -> 0.291575}}, {1.51666, {zQ -> 
     0.38433}}, {0.936696, {zQ -> 0.476557}}, {0.636637, {zQ -> 
     0.568672}}, {0.461284, {zQ -> 0.660829}}, {0.349854, {zQ -> 
     0.753076}}, {0.274588, {zQ -> 0.845422}}}}

In a case of n=1 we can plot

Plot[Evaluate[
  Table[torootsig[.1, x, y, .1], {y, .093, .0999, .0011}]], {x, -.01, 
  0.01}, PlotRange -> {-.05, .05}]

Figure 1

There are no roots in a range of 0.093<=xQ<=0.098 it is why we have messages from FindRoot, NIntegrate and FindMinimum. To solve this problem we can use

cp = ContourPlot[
  torootsig[.1, x, y, .1] == 0, {x, -.1, 0.0}, {y, .093, .0999}]

Figure 2

To retrieve data from cp we evaluate

data = cp[[1, 1]][[1]]

Out[]= {{0., 0.0990922}, {-0.000597747, 0.0991195}, {-0.0014888, 
  0.0991607}, {-0.00279017, 0.0992146}, {-0.00357143, 
  0.0992546}, {-0.00417006, 0.0992839}, {-0.00487564, 
  0.0993172}, {-0.00535714, 0.0993408}, {-0.00590494, 
  0.0993693}, {-0.00659417, 0.0994071}, {-0.00690503, 
  0.0994236}, {-0.00714286, 0.0994372}, {-0.00786928, 
  0.0994802}, {-0.00863636, 0.0995304}, {-0.00878632, 
  0.0995402}, {-0.00892857, 0.0995513}, {-0.00946648, 
  0.099592}, {-0.00963374, 0.0996049}, {-0.00982143, 
  0.0996235}, {-0.0101356, 0.0996536}, {-0.0103567, 
  0.0996782}, {-0.0105918, 0.0997152}, {-0.0107143, 
  0.0997582}, {-0.01076, 0.0997736}, {-0.0107628, 
  0.0997768}, {-0.010754, 0.0997795}, {-0.0107143, 
  0.0997873}, {-0.0105213, 0.0998384}, {-0.0102057, 
  0.0998649}, {-0.00982143, 0.0998815}, {-0.00958534, 0.0999}}

These data can be used to minimize SQ in a case of n=1,2.

$\endgroup$
8
  • $\begingroup$ Nice analysis of the problem, however, I'm already aware of the non-existence of the root there, that is why I'm asking if there is any way to procedurally change the range of FindMinimum to avoid the problem and generate a complete list of the minimum point in a single Table. Your approach is also similar to my original approach to this problem, kind of manually searching for where FindRoot will not fail. However, I'm thinking of a general case so that if I go even below zh<0.1, I will have a procedure that bypasses this problem so that I don't to do it for every n. $\endgroup$
    – mathemania
    Aug 31, 2021 at 7:37
  • $\begingroup$ There is no any general equation for your problem it is why I try to keep your code. May be we can compute your problem in one line. Can you explain what do you try to calculate? $\endgroup$ Aug 31, 2021 at 8:53
  • $\begingroup$ As in my post above, I'm calculating the point zQ where SQ is a minimum. I want to do that for some zh. Now I can do it one by one for every zh but as you see zQmin of the range {zQ, zQmin, zQmax} will need to be increased for lower zh. I already have the analysis of the problem why it does not work from the link I pasted in my post and your answer. But now, what if I want to go lower than zh=0.1, what if I will compute for {zh,0.01,0.09,0.01}? So I'm thinking maybe there is a general scheme for this problem based on the pattern of the solution above. $\endgroup$
    – mathemania
    Aug 31, 2021 at 9:31
  • $\begingroup$ I used your first plot to study for example zh=0.03, it seems like it does not cross the x-axis even for 0.9999999 zh. So maybe there is a zh where torootsig does not have a root already. Plot[torootsig[0.1, x, 0.9999999 0.03, 0.03], {x, -0.005, 0.005}, PlotRange -> {-0.05, 0.05}, PlotStyle -> {Blue, Thickness[0.005]}, ImageSize -> Large]. The question is how to find this zh where torootsig fails to have a root. $\endgroup$
    – mathemania
    Aug 31, 2021 at 10:56
  • $\begingroup$ @mathemania Your code is based on definition sig[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := x /. FindRoot[torootsig[b, x, zQ, zh], {x, -2 zh, 0}, MaxIterations -> 1000]; while there are no roots for some parameters. Therefore, first, we need to separate regions with and without roots. Second, we need to define what to do in the case when there are no roots. $\endgroup$ Sep 1, 2021 at 1:12

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.