0
$\begingroup$

I have a function SQ[b, zQ, zh] where I can find at which zQ it is a minimum, b=0.1 and zh can be specified. Similarly, using the zQ I found, I can use it to find the value of sig[b, zQ, zh] with the same parameters b and zh. Again, the value of zQ is found corresponding to the zh specified.

Next, I have a function SvQ[zm, zQ, sig, zh] where we use the determined values of zQ and sig above, zh is the same specified value above, so that SvQ[zm] varies as a function of zm. An example plot of SvQ[zm] for zh=0.5 (left) and zh=0.6 (right) are shown below.

Image

You can see that the minimum value of SvQ[zm] for zh=0.5 is a negative value, while for zh=0.6 it is a positive value.

I want to find the zh so that SvQ[zm] has a minimum value of 0, or at least essentially 0, like 10^-10 or maybe even higher.

The issue now is that I do this manually, i.e. I choose zh to vary from 0.5600 to 0.5700 (the minimum of SvQ is negative for 0.5600 and positive for 0.5700) by steps of 10, but this will only tell me the value of zh with only four decimal places. You can see in FindMinimum[SvQ] below that it gave me a result that zh=0.5650 is negative and zh=0.5660 is positive, I can continue increasing the digits manually but very inefficient.

Is there a better way of determining at which zh the minimum of SvQ is equal to 0?

ClearAll["Global`*"]
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 -> 100]]
sig[b_?NumericQ, zQ_?NumericQ, zh_?NumericQ] := sig /. FindRoot[torootsig[b, sig, zQ, zh] == 0, {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 ]

test1 = ParallelTable[FindMinimum[SQ[0.1, zQ, zh/10^4], {zQ, (9950/10^4) (zh/10^4), (9300/10^4) (zh/10^4), (9990/10^4) (zh/10^4)}, AccuracyGoal -> ag/2, PrecisionGoal -> pg/2, WorkingPrecision -> wp, MaxIterations -> 100], {zh, 5600, 5700, 10}] // AbsoluteTiming
test1List = Table[zQ /. Extract[Table[test1[[2, i, 2]], {i, 1, 11}], {i, 1}], {i, 1, 11}]
test2List = MapThread[sig[1/10, #1, (#2/10^4)] &, {test1List, Table[i, {i, 5600, 5700, 10}]}]

SvQ1[z_?NumericQ, zm_?NumericQ, zQ_?NumericQ, sig_?NumericQ, zh_?NumericQ] := Module[{zr, zmr, zQr, sigr, zhr}, {zr, zmr, zQr, sigr, zhr} = Rationalize[{z, zm, zQ, sig, zh}, 0]; zmr^d/(4 zr^d Sqrt[f[zr, zhr] zmr^(2 d) - f[zmr, zhr] zr^(2 d)]) - zQr^d/(4 zr^d Sqrt[f[zr, zhr] (zQr^(2 d) - (f[zQr, zhr] zr^(2 d))/(f[zQr, zhr] + sigr^2))])]
SvQ2[z_?NumericQ, zm_?NumericQ, zh_?NumericQ] := Module[{zr, zmr, zhr}, {zr, zmr, zhr} = Rationalize[{z, zm, zh}, 0]; zmr^d/(4 zr^d Sqrt[f[zr, zhr] zmr^(2 d) - f[zmr, zhr] zr^(2 d)])]
SvQ[zm_?NumericQ, zQ_?NumericQ, sig_?NumericQ, zh_?NumericQ] := Module[{zmr, zQr, sigr, zhr}, {zmr, zQr, sigr, zhr} = Rationalize[{zm, zQ, sig, zh}, 0]; NIntegrate[SvQ1[z, zmr, zQr, sigr, zhr], {z, 0, zQr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 100] + 1/(4 zmr^(d - 1)) - 1/(4 zQr^(d - 1)) + NIntegrate[SvQ2[z, zmr, zhr], {z, zQr, zmr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 100]]

MapThread[FindMinimum[SvQ[zm, #1, #2, #3/10^4], {zm, 1.03725505 (#3/10^4), 1.31605 (#3/10^4)}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 100] &, {test1List, test2List, Table[i, {i, 5600, 5700, 10}]}] // AbsoluteTiming

{402.8128997,{{-0.0007699902743,{zm->0.6427891763}},{-0.0006186271438,{zm->0.6439370142}},{-0.0004690419362,{zm->0.645084852}},{-0.0003212141805,{zm->0.6462326898}},{-0.0001751236627,{zm->0.6473805276}},{-0.00003075042224,{zm->0.6485283654}},{0.0001119252514,{zm->0.6496762032}},{0.0002529228224,{zm->0.6508240411}},{0.0003922615119,{zm->0.6519718789}},{0.0005299603013,{zm->0.6531197167}},{0.0006660379357,{zm->0.6542675545}}}}

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.