I need to calculate Sv[zm,zs,zh]
for a given range of zm
and zs
with fixed zh
. I have written a code where zm
is a function of zs
given by zmzs[zs, zh, tb, zmmin, zmmax]
so that Sv[zm(zs),zs,zh] = Sv[zmzs[zs, zh, tb, zmmin, zmmax],zs,zh]
and now I only need to give the range of zs
only. To add, zh
and tb
are given constants, zmmin
and zmmax
are the min and max of the range in FindRoot
in zmzs[zs, zh, tb, zmmin, zmmax]
.
d = 3;
ag = 10;
pg = 10;
wp = 20;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
tzsint[z_?NumericQ, zs_?NumericQ, zh_?NumericQ] := Module[{zr, zsr, zhr}, {zr, zsr, zhr} = Rationalize[{z, zs, zh}, 0]; zr^d/Sqrt[f[zr, zhr] (zsr^(2 d) - zr^(2 d))]]
tzs[zs_?NumericQ, zh_?NumericQ] := Module[{zsr, zhr}, {zsr, zhr} = Rationalize[{zs, zh}, 0]; NIntegrate[tzsint[z, zsr, zhr], {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20]]
tzmint[z_?NumericQ, zm_?NumericQ, zh_?NumericQ] := Module[{zr, zmr, zhr}, {zr, zmr, zhr} = Rationalize[{z, zm, zh}, 0]; -1/(f[zr, zhr] Sqrt[1 - (zmr^(2 d) f[zr, zhr])/(z^(2 d) f[zmr, zhr])])]
tzm[zm_?NumericQ, zh_?NumericQ] := Module[{zmr, zhr}, {zmr, zhr} = Rationalize[{zm, zh}, 0]; NIntegrate[tzmint[z, zmr, zhr], {z, 0, zhr, zmr}, Method -> PrincipalValue, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20]]
Svint1[z_?NumericQ, zm_?NumericQ, zh_?NumericQ] := Module[{zr, zmr, zhr}, {zr, zmr, zhr} = Rationalize[{z, zm, zh}, 0]; zmr^d/(zr^d Sqrt[f[zr, zhr] zmr^(2 d) - f[zmr, zhr] zr^(2 d)])]
Svint2[z_?NumericQ, zm_?NumericQ, zs_?NumericQ, zh_?NumericQ] := Module[{zr, zmr, zsr, zhr}, {zr, zmr, zsr, zhr} = Rationalize[{z, zm, zs, zh}, 0]; (((zsr^(2 d) zmr^(2 d) zr)/zhr^(d + 1)) + (zr^d (zsr^(2 d) f[zmr, zhr] - zmr^(2 d))))/((zmr^d Sqrt[zsr^(2 d) - zr^(2 d) ] + zsr^d Sqrt[zmr^(2 d) f[zr, zhr] - zr^(2 d) f[zmr, zhr]]) Sqrt[(zmr^(2 d) f[zr, zhr] - zr^(2 d) f[zmr, zhr]) (zsr^(2 d) - zr^(2 d))])]
Sv[zm_?NumericQ, zs_?NumericQ, zh_?NumericQ] := Module[{zmr, zsr, zhr}, {zmr, zsr, zhr} = Rationalize[{zm, zs, zh}, 0]; NIntegrate[Svint1[z, zmr, zhr], {z, zsr, zmr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20] + NIntegrate[Svint2[z, zmr, zsr, zhr], {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20]]
tzmzs[zm_?NumericQ, zs_?NumericQ, zh_?NumericQ, tb_?NumericQ] := Module[{zmr, zsr, zhr, tbr}, {zmr, zsr, zhr, tbr} = Rationalize[{zm, zs, zh, tb}, 0]; NIntegrate[tzmint[z, zmr, zhr], {z, 0, zhr, zmr}, Method -> PrincipalValue, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20] - NIntegrate[tzsint[z, zsr, zhr], {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20] - tbr]
zmzs[zs_?NumericQ, zh_?NumericQ, tb_?NumericQ, zmmin_?NumericQ, zmmax_?NumericQ] := zm /. FindRoot[tzmzs[zm, zs, zh, tb] == 0, {zm, zmmin, zmmax}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
zszm[zm_?NumericQ, zh_?NumericQ, tb_?NumericQ, zsmin_?NumericQ, zsmax_?NumericQ] := zs /. FindRoot[tzmzs[zm, zs, zh, tb] == 0, {zs, zsmin, zsmax}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20]
The issue is calculating zmzs[zs, zh, tb, zmmin, zmmax]
takes a lot of time, for example I evaluated everything at zs=9.9952388099
. zmzs[zs, zh, tb, zmmin, zmmax]
takes around 13 secs to calculate, Sv[zmzs[zs, zh, tb, zmmin, zmmax],zs,zh]
takes 15 secs to calculate. Using the value of zm=13.127957300691564691
calculated from zmzs[zs, zh, tb, zmmin, zmmax]
and plugging it directly to Sv[zm,zs,zh]
takes only around 2 secs. This means that Sv
itself calculates fairly quick, it is really zmzs[zs, zh, tb, zmmin, zmmax]
that takes quite some time.
In the end I want to plot Sv
for [zs,1.5,9.9952388097]
which will definitely take a lot of time, in fact I tried plotting it and it has already been 4 hrs but so far no result yet!!!
zmzs[9.9952388099, 10, 0.100000100000003174, 10.00099982369714445, 13.1279573]//AbsoluteTiming
{13.1723073, 13.127957300691564691}
Sv[zmzs[9.9952388099, 10, 0.100000100000003174, 10.00099982369714445, 13.1279573], 9.9952388099, 10]//Chop//AbsoluteTiming
{15.0956025, 0.010820232314756221363}
Sv[13.127957300691564691, 9.9952388099, 10]//Chop//AbsoluteTiming
{1.9989049, 0.010820232314756221363}
Is there any way to speed up the code in order to plot out Sv
in a shorter amount of time?
ListPlot[..,Joined->True]
and calculate just enough points for that to be able to see the behavior without needing to have Plot calculate every single pixel. $\endgroup$WorkingPrecision
above the standard value. It would be nice to know why. $\endgroup$