Wygładź wynik FindRoot

Question

Mam zestaw kodu, dla którego wymaga to znalezienia odpowiedniego cdla każdego a(chociaż podam wartość apóźniej) i zużycia ograniczenia, toroot[a,c,z]a następnie podstawienia z cpowrotem do końcowego wyrażenia functionS[a,z]. Mam też inną funkcję, dla której jest zmiana zmiennej functionSR[l,z]gdzie 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];

Mój problem jest, gdy próbuję znaleźć minimum functionS[a,z]i functionSR[l,z]dla niektórych ai l, powiedzmy a=1, a l=1to daje mi błąd. Myślę, że jest to związane z zachowaniem, ckiedy a=1lub 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}]}

Dla a = 0,1 wykres jest znacznie bardziej płynny

Dla a = 1 wykres zawiera więcej wypukłości

Czy mój kod jest źle napisany do wyodrębnienia c? Czy są jakieś zmiany, które można wprowadzić? Czytałem gdzieś, że Reducemożna go również użyć zamiast, FindRootale wciąż to rozumiem. Czy użycie LocalAdaptivejako metody jest NIntegratetutaj odpowiednie?

AKTUALIZACJA: Zwróć uwagę na literówkę, poprawiłem ją. Na działkach wcześniej pisałem c=0.1i c=1ale powinno być a=0.1i a=1.

Wyrażenia mojego problemu są podane przez:

$$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}

gdzie $b=\frac{z_s}{z_h}$ i zanotuj to $c=c(z_s)$( c=c[z]) chociaż w kodzie c=c[a,z],$c$ powinno zależeć tylko od $z_s$( z) od$a$ zostanie określony na końcu.

Być może jest też lepszy sposób na znalezienie projektu $c$. Właściwie mogę mieć inne ograniczenie, w którym$\frac{dS}{dz_s} = 0$ (to dlatego, że w końcu muszę zminimalizować $S$ z szacunkiem do $z_s$) i być może pochodną $\eqref{1}$ z szacunkiem do $z_s$, aby można je było znaleźć $c$?

Answer 1

Źródło NIntegratekomunikatów o błędach można zobaczyć na podstawie współczynnika całkowania,, x^d/Sqrt[1-c x^d z^d]z toroot. Dla c > z^-3The podcałkowa jest w liczbie pojedynczej na jakimś punkcie w domenie {x, 0, 1}. Co więcej, gdyby NIntegratemożna było całkować poprzez osobliwość (a z pomocą może), wynikiem byłaby liczba zespolona, ​​co (przypuszczalnie) jest niepożądane. Aby kontynuować, zmień zmienną całkowania na xd = x^(d+1)i zastosuj odpowiednią Methodz tego miejsca .

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

Ponadto przedefiniuj, czaby używać siecznego Methodi ogranicz wyszukiwanie cmiędzy 0i z^-3.

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

(Początkowe domysły .5 z^-3i .6 z^-3, zostały wybrane nieco arbitralnie.) Z tą definicją czzwraca poprawną wartość c, jeśli istnieje, i z^-3wraz z FindRoot::regedkomunikatem o błędzie w przeciwnym razie. Dzięki tym definicjom oba wykresy w pytaniu można uzyskać poprawnie w następujący sposób. dla a = 1,

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

Checkzapobiega wykreślaniu krótkiego zasięgu wokół c = 1.42, gdzie nie ma rozwiązania, chociaż nie eliminuje odpowiednich komunikatów o błędach. Drugi wątek a = .1to

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