0
$\begingroup$

I have a code where the first line of the body consists of an algebraic equation involving NIntegrate and the parameter a, for a given a the solution will be stored as a function toroot[c_?NumericQ, z_], c will then be solved as a function of z and substituted to the expression functionS.

This is all good and I initially placed a value for a=0.2, but

(1) What strategy should I employ so that I can write a functionS where I can plug in several values for a so that in the end I can plot functionS vs. z_s for different a?

(2) I want to find the minimum of functionS using FindMinimum for different a, that is, is there a way to do this in a single code or maybe produce a list of the minimum values. Also, I want to plot the minimum points as a function of a.

I believe I will use Table but I'm not quite sure how to write the proper efficient code.

d = 3;
zh = 1.5;
a = 0.2;
toroot[c_?NumericQ, z_] := a - NIntegrate[(c z^(d + 1) x^d)/((1 - ((z x)/zh)^(d + 1)) (1 - c^2 (z x)^(2 d)))^(1/2), {x, 0, 1}, MaxRecursion -> 5, PrecisionGoal -> 4, Method -> "DuffyCoordinates"]
cz[z_?NumericQ] := c /. FindRoot[toroot[c, z], {c, 0.009, 0.0000001, 100}, WorkingPrecision -> 3]
intS[z_?NumericQ] := NIntegrate[With[{b = z/1.5}, (((-1)/(d - 1)) cz[z]^2 z^(2 d)) x^d ((1 - (b x)^(d + 1))/(1 - cz[z]^2 (z x)^(2 d)))^(1/2) - ((b^(d + 1) (d + 1))/(2 (d - 1))) x ((1 - cz[z]^2 (z x)^(2 d))/(1 - (b x)^(d + 1)))^(1/2) + (b^(d + 1) x)/((1 - (b x)^(d + 1)) (1 - cz[z]^2 (z x)^(2 d)))^(1/2)], {x, 0, 1}, MaxRecursion -> 5, PrecisionGoal -> 4, Method -> "DuffyCoordinates"]
functionS = ((-((1 - cz[z]^2 z^(2 d)) (1 - (z/1.5)^(d + 1)))^(1/2)/(d - 1)) + intS[z] + 1)/(4 z^(d - 1));
function = Log[10, functionS];

Plot[function, {z, 1.3, 1.5}, PlotPoints -> 3, AxesLabel -> {"z_s", "Log S"}, AxesStyle -> Directive[Black, 18], PlotStyle -> {Thick, Blue}, PlotRange -> Full, ImageSize -> Large] // Quiet // AbsoluteTiming

FindMinimum[functionS, {z, 1.46}] // Quiet // AbsoluteTiming

Pics

This plot is only functionS for a=0.2. What I want to get is a plot of functionS say for a=0.1,0.2,0.3,0.4,0.5.

UPDATE: Question (1) is already resolved (the reason I cannot do the plots is because I defined toroot[a_?NumericQ,c_?NumericQ, z_] instead of toroot[a_,c_?NumericQ, z_] as NIntegrate would require a numerical value), BUT (2) is not yet, can anyone help me on this?

$\endgroup$
7
  • $\begingroup$ Just define everything as a function of a? $\endgroup$ – Natas Sep 30 '20 at 16:10
  • $\begingroup$ @Natas I have tried that but failed but now I already know why your suggestion did not work, because I need to write toroot[a_,c_?NumericQ, z_] instead of toroot[a_?NumericQ,c_?NumericQ, z_]. However, my question about FindMinimum is still unanswered I want to write a code that gives a minimum of functionS for different a, also I want to plot the minimum points as a function of a. $\endgroup$ – mathemania Sep 30 '20 at 16:38
  • $\begingroup$ To find the minima: minima = Table[First@FindMinima[functionS[a], {z, 1.3,1.5}], {a, 0.1, 0.5, 0.1}]; To plot them as a function of a: ListPlot[minima, DataRange -> {0.1, 0.5}]. $\endgroup$ – Natas Sep 30 '20 at 19:39
  • $\begingroup$ @Natas Your suggestion worked! However, do you know if my code can be made any faster? Since it takes a long time just to plot say, 5 points. I have used DuffyCoordinates in NIntegrate to avoid any singularities and maybe speed up the calculation. I have also introduced WorkingPrecision and PrecisionGoal, also I added MaxRecursion but I am not quite sure about its real use, any thoughts? $\endgroup$ – mathemania Oct 1 '20 at 15:15
  • $\begingroup$ Unfortunately I have very little experience with numerics in Mathematica. Certainly decreasing the options you mentioned can increase the computation speed. However, of course at a loss of precision. I had a look at the documentation of "DuffyCoordinates" and don't think you should use it. Firstly, it's for multidimensional integration and, secondly, you also need to specify "Corners" (relevant documentation). $\endgroup$ – Natas Oct 2 '20 at 8:44

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.