Mathematica Asked by mathemania on December 6, 2020

I have a set of code for which it involves finding the corresponding `c`

for each `a`

(although I will give a value of `a`

later) and `z`

using the constraint `toroot[a,c,z]`

and then substituting `c`

back into the final expression `functionS[a,z]`

. I also have another function for which there is a change of variable `functionSR[l,z]`

where `a->l-0.01`

.

```
d = 3;
zh = 1.5;
toroot[a_, 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 -> "LocalAdaptive"]
cz[a_?NumericQ, z_?NumericQ] := c /. FindRoot[toroot[a, c, z], {c, 0.0009, 0.0000001, 10000}, WorkingPrecision -> 5]
intS[a_?NumericQ, z_?NumericQ] := NIntegrate[With[{b = z/zh}, (((-1)/(d - 1)) cz[a, z]^2 z^(2 d)) x^d ((1 - (b x)^(d + 1))/(1 - cz[a, z]^2 (z x)^(2 d)))^(1/2) - ((b^(d + 1) (d + 1))/(2 (d - 1))) x ((1 - cz[a, z]^2 (z x)^(2 d))/(1 - (b x)^(d + 1)))^(1/2) + (b^(d + 1) x)/((1 - (b x)^(d + 1)) (1 - cz[a, z]^2 (z x)^(2 d)))^(1/2)], {x, 0, 1}, MaxRecursion -> 5, PrecisionGoal -> 4, Method -> "LocalAdaptive"]
functionS[a_, z_] = ((-((1 - cz[a, z]^2 z^(2 d)) (1 - (z/zh)^(d + 1)))^(1/2)/(d - 1)) + intS[a, z] + 1)/(z^(d - 1));
functionSR[l_, z_] = Replace[functionS[a, z], a -> (l - 0.01), Infinity];
```

My problem is when I try to find the minimum of `functionS[a,z]`

and `functionSR[l,z]`

for some `a`

and `l`

, say `a=1`

and `l=1`

, it gives me an error. I think it is connected to the behavior of `c`

when `a=1`

or `l=1`

.

```
In[23]:= FindMinimum[functionS[1, z], {z, 1.2, 1.5}] //
Quiet // AbsoluteTiming
FindMinimum[functionSR[1, z], {z, 1.2, 1.5}] // Quiet // AbsoluteTiming
During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000944548 -0.00149313 I and 0.0006178735732839699` for the integral and error estimates.
During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000944548 -0.00149313 I and 0.0006178735732839699` for the integral and error estimates.
During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000949747 -0.00149122 I and 0.000620731102746343` for the integral and error estimates.
During evaluation of In[23]:= General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation.
During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.
During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.
During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.
During evaluation of In[23]:= General::stop: Further output of FindRoot::reged will be suppressed during this calculation.
During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.436961 -1.38189 I is not a real number at {z} = {1.2}.
During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.436961 -1.38189 I is not a real number at {z} = {1.2}.
Out[23]= {0.760891, FindMinimum[functionS[1, z], {z, 1.2, 1.5}]}
During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286247 -0.0000971587 I and 0.0005426332486649041` for the integral and error estimates.
During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286247 -0.0000971587 I and 0.0005426332486649041` for the integral and error estimates.
During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286812 -0.0000961916 I and 0.0005442259497809905` for the integral and error estimates.
During evaluation of In[23]:= General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation.
During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.
During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.
During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.
During evaluation of In[23]:= General::stop: Further output of FindRoot::reged will be suppressed during this calculation.
During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.439434 -1.36539 I is not a real number at {z} = {1.2}.
During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.439434 -1.36539 I is not a real number at {z} = {1.2}.
Out[24]= {0.771827, FindMinimum[functionSR[1, z], {z, 1.2, 1.5}]}
```

For a=0.1, the plot is much more smooth

For a=1, the plot contains more bumps

Is my code badly written to extract `c`

? Are there any changes that can be done? I have read somewhere that `Reduce`

can also be used instead of `FindRoot`

but I am still figuring it out. Also, is using `LocalAdaptive`

as a method for `NIntegrate`

suitable here?

UPDATE:

Please note of the typo, I have corrected it. In the plots before, I wrote `c=0.1`

and `c=1`

but it should be `a=0.1`

and `a=1`

.

The expressions of my problem is given by,

$$a = c z_s^{d+1}int_0^1 dx frac{x^d}{sqrt{(1-(z_s/z_h)^{d+1} x^{d+1})(1-c^2 z_s^{2d} x^{2d})}} tag{1}label{1}$$

begin{align}

S &= frac{1}{4 z_s^{d-1}}Bigg(1 -frac{sqrt{(1-c^2 z_s^{2d})(1-b^{d+1})}}{d-1} – frac{1}{d-1} c^2 z_s^{2d} int^1_0 dx x^d sqrt{frac{(1-(b x)^{d+1})}{(1-c^2(z_s x)^{2d})}}\

& -frac{b^{d+1}(d+1)}{2(d-1)} int^1_0 dx x sqrt{frac{(1-c^2(z_s x)^{2d})}{(1-(b x)^{d+1})}}\

& + b^{d+1}int^1_0 dx frac{x}{sqrt{(1-(b x)^{d+1})(1-c^2(z_s x)^{2d})}}Bigg) tag{2}label{2}

end{align}

where $b=frac{z_s}{z_h}$ and note that $c=c(z_s)$ (`c=c[z]`

) although in the code `c=c[a,z]`

, $c$ should only depend on $z_s$ (`z`

) since $a$ will be specified in the end.

Also, maybe there is a better way to design finding $c$. Actually, I can have another constraint where $frac{dS}{dz_s} = 0$ (that is because in the end I need to minimize $S$ with respect to $z_s$) and maybe the derivative of $eqref{1}$ with respect to $z_s$, so that these can be used to find $c$?

The source of the `NIntegrate`

error messages can be seen from a factor of the integrand, `x^d/Sqrt[1-c x^d z^d]`

, of `toroot`

. For `c > z^-3`

, the integrand is singular for some point in the domain, `{x, 0, 1}`

. Moreover, if `NIntegrate`

could integrate through the singularity (and, with help, it can), the result would be a complex number, which (presumably) is undesirable. To proceed, change the variable of integration to `xd = x^(d+1)`

and apply the appropriate `Method`

from here.

```
toroot[a_, c_?NumericQ, z_] := a - NIntegrate[((1 - xd (z /zh)^(d + 1))
(1 - c^2 xd^(2 d/(d + 1)) z^(2 d)))^(-1/2), {xd, 0, 1}, Method -> {"GlobalAdaptive",
"SingularityHandler" -> "DoubleExponential"}] (c z^(d + 1))/4
```

In addition, redefine `cz`

to use the secant `Method`

and bound the search for `c`

to between `0`

and `z^-3`

.

```
cz[a_?NumericQ, z_?NumericQ] := c /.
FindRoot[toroot[a, c, z], {c, .5 z^-3, .6 z^-3/2, 0, z^-3}]
```

(The initial guesses, `.5 z^-3`

and `.6 z^-3`

, were chosen somewhat arbitrarily.) With this definition, `cz`

returns the correct value of `c`

, if it exists, and `z^-3`

along with the `FindRoot::reged`

error message otherwise. With these definitions, the two plots in the question can be obtained correctly as follows. For `a = 1`

,

```
Plot[Check[cz[1, z], Null], {z, 1.42, zh}, AxesLabel -> {z, c},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
```

`Check`

prevents plotting of a short range around `c = 1.42`

, where no solution exists, although it does not eliminate the corresponding error messages. The second plot, for `a = .1`

, is

```
LogPlot[Check[cz[.1, z], Null], {z, .2, zh}, AxesLabel -> {z, c},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
```

Correct answer by bbgodfrey on December 6, 2020

Get help from others!

Recent Answers

- Joshua Engel on Why fry rice before boiling?
- Peter Machado on Why fry rice before boiling?
- Jon Church on Why fry rice before boiling?
- haakon.io on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?

Recent Questions

- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP