1
$\begingroup$

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ How fast is it if you get rid of all the Rationalize and all the non-machine-precision calculations? Then all the calculations should run with your CPU hardware floating point and perhaps be substantially faster. Then how fast is it if instead of Plot your use 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$
    – Bill
    Sep 5, 2021 at 13:23
  • $\begingroup$ Can you add equations you are trying compute. It seems you do multiple numerical integrations using WorkingPrecision above the standard value. It would be nice to know why. $\endgroup$
    – yarchik
    Sep 5, 2021 at 15:03

1 Answer 1

Reset to default
2
$\begingroup$ 
Clear["Global`*"]

d = 3;
ag = 10;
pg = 10;
wp = 20;

SetOptions[NIntegrate,
  AccuracyGoal -> ag,
  PrecisionGoal -> pg,
  WorkingPrecision -> wp,
  MaxRecursion -> 20];

SetOptions[FindRoot,
  AccuracyGoal -> ag,
  PrecisionGoal -> pg,
  WorkingPrecision -> wp];

f[z_, zh_] := 1 - (z/zh)^(d + 1);

It is only necessary to restrict a function to numeric arguments if the function makes use of a numeric technique.

tzsint[z_, zs_, zh_] = z^d/Sqrt[f[z, zh] (zs^(2 d) - z^(2 d))];

tzs[zs_?NumericQ, zh_?NumericQ] :=
 Module[{zsr, zhr},
  {zsr, zhr} = Rationalize[{zs, zh}, 0];
  NIntegrate[tzsint[z, zsr, zhr], {z, 0, zsr}]]

Complicated calculations that do not use numeric techniques should be done once (Set vice SetDelayed or use Evaluate on RHS of SetDelayed) and use Simplify to help improve efficiency (EDIT) by reducing the number of individual operations.

tzmint[z_, zm_, zh_] =
  -1/(f[z, zh] Sqrt[1 - (zm^(2 d) f[z, zh])/(z^(2 d) f[zm, zh])]) //
   Simplify;

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]]

Svint1[z_, zm_, zh_] =
  zm^d/(z^d Sqrt[f[z, zh] zm^(2 d) - f[zm, zh] z^(2 d)]) //
   Simplify;

Svint2[z_, zm_, zs_, zh_] =
  (((zs^(2 d) zm^(2 d) z)/zh^(d + 1)) + (z^d (zs^(2 d) f[zm, zh] - zm^(2 d))))/
    ((zm^d Sqrt[zs^(2 d) - z^(2 d)] + 
        zs^d Sqrt[
          zm^(2 d) f[z, zh] - z^(2 d) f[zm, zh]]) Sqrt[(zm^(2 d) f[z, zh] - 
          z^(2 d) f[zm, zh]) (zs^(2 d) - z^(2 d))]) //
   Simplify;

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}] +
   NIntegrate[Svint2[z, zmr, zsr, zhr], {z, 0, zsr}]]

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] -
   NIntegrate[tzsint[z, zsr, zhr], {z, 0, zsr}] - tbr]

zmzs[zs_?NumericQ, zh_?NumericQ, tb_?NumericQ, zmmin_?NumericQ, 
  zmmax_?NumericQ] :=
 zm /. FindRoot[tzmzs[zm, zs, zh, tb] == 0,
   {zm, zmmin, zmmax}]

zszm[zm_?NumericQ, zh_?NumericQ, tb_?NumericQ, zsmin_?NumericQ, 
  zsmax_?NumericQ] :=
 zs /. FindRoot[tzmzs[zm, zs, zh, tb] == 0,
   {zs, zsmin, zsmax}]

Checking the timing

zmzs[9.9952388099, 10, 0.100000100000003174, 10.00099982369714445, 
  13.1279573] //
 AbsoluteTiming

(* {0.272723, 13.127957300691564608} *)

Sv[zmzs[9.9952388099, 10, 0.100000100000003174, 10.00099982369714445, 
    13.1279573], 9.9952388099, 10] // Chop // AbsoluteTiming

(* {0.276597, 0.0108202323147562170881} *)

Sv[13.127957300691564691, 9.9952388099, 10] // Chop // AbsoluteTiming

(* {0.015935, 0.010820232314756221330} *)

You did not indicate the range of interest for zs for the plot

Plot[
   Sv[zs, 1.5, 9.9952388097],
   {zs, 3/2, 15},
   PlotRange -> All,
   PlotPoints -> 50,
   MaxRecursion -> 5,
   ImageSize -> 400] //
  AbsoluteTiming //
 Column

enter image description here

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ Re "use Simplify to help preserve precision by reducing the number of individual operations": Error propagation is not a concern of Simplify, which is mainly concerned with reducing the number of leaves of the expression. I don't know how often — I would guess not very — but it sometimes makes the error much worse. But overall, +1. Just didn't want others to think Simplify can be relied on in this way. $\endgroup$
    – Michael E2
    Sep 5, 2021 at 16:07
  • $\begingroup$ Example: mathematica.stackexchange.com/a/236767 -- as I said, I expect this is a rare case, but Simplify makes no attempt to avoid underflow or subtractive cancellation. $\endgroup$
    – Michael E2
    Sep 5, 2021 at 16:17
  • $\begingroup$ @MichaelE2 - Thanks. Corrected. $\endgroup$
    – Bob Hanlon
    Sep 5, 2021 at 16:23

Your Answer

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

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