Uniforma il risultato di FindRoot
Ho un set di codice per il quale si tratta di trovare il corrispondente cper ciascuno a(anche se darò un valore di apiù tardi) e di zutilizzare il vincolo toroot[a,c,z]e quindi di sostituirlo cnell'espressione finale functionS[a,z]. Ho anche un'altra funzione per la quale c'è un cambio di variabile functionSR[l,z]dove 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];
Il mio problema è quando cerco di trovare il minimo di functionS[a,z]e functionSR[l,z]per alcuni ae l, diciamo a=1e l=1, mi dà un errore. Penso che sia collegato al comportamento di cquando a=1o 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}]}
Per a = 0,1, il grafico è molto più fluido
Per a = 1, il grafico contiene più dossi
Il mio codice è scritto male per estrarre c? Ci sono cambiamenti che possono essere fatti? Ho letto da qualche parte che Reducepuò essere utilizzato anche al posto di, FindRootma lo sto ancora cercando di capire. Inoltre, sta usando LocalAdaptivecome metodo NIntegrateadatto qui?
AGGIORNAMENTO: Si prega di notare l'errore di battitura, l'ho corretto. Nelle trame prima, ho scritto c=0.1e c=1ma dovrebbe essere a=0.1e a=1.
L'espressione del mio problema è data da,
$$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}
dove $b=\frac{z_s}{z_h}$ e nota quello $c=c(z_s)$( c=c[z]) sebbene nel codice c=c[a,z],$c$ dovrebbe dipendere solo da $z_s$( z) da$a$ sarà specificato alla fine.
Inoltre, forse c'è un modo migliore per trovare il design $c$. In realtà, posso avere un altro vincolo dove$\frac{dS}{dz_s} = 0$ (questo perché alla fine ho bisogno di ridurre al minimo $S$ riguardo a $z_s$) e forse il derivato di $\eqref{1}$ riguardo a $z_s$, in modo che possano essere utilizzati per trovare $c$?
Risposte
La fonte dei NIntegratemessaggi di errore può essere vista da un fattore dell'integrando,, x^d/Sqrt[1-c x^d z^d]di toroot. Per c > z^-3, l'integrando è singolare per un certo punto del dominio, {x, 0, 1}. Inoltre, se NIntegratepotesse integrarsi attraverso la singolarità (e, con l'aiuto, può), il risultato sarebbe un numero complesso, che (presumibilmente) è indesiderabile. Per procedere, modificare la variabile di integrazione xd = x^(d+1)e applicare l'appropriata Methodda qui .
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
Inoltre, ridefinisci czper utilizzare la secante Methode lega la ricerca ca tra 0e z^-3.
cz[a_?NumericQ, z_?NumericQ] := c /.
FindRoot[toroot[a, c, z], {c, .5 z^-3, .6 z^-3/2, 0, z^-3}]
(Le ipotesi iniziali, .5 z^-3e .6 z^-3, sono state scelte in modo piuttosto arbitrario.) Con questa definizione, czrestituisce il valore corretto di c, se esiste, e z^-3insieme al FindRoot::regedmessaggio di errore altrimenti. Con queste definizioni, i due grafici nella domanda possono essere ottenuti correttamente come segue. per a = 1,
Plot[Check[cz[1, z], Null], {z, 1.42, zh}, AxesLabel -> {z, c},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
Checkimpedisce la stampa di un breve raggio intorno c = 1.42, dove non esiste soluzione, sebbene non elimini i messaggi di errore corrispondenti. La seconda trama, per a = .1, è
LogPlot[Check[cz[.1, z], Null], {z, .2, zh}, AxesLabel -> {z, c},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
