Сгладить результат FindRoot

Question

У меня есть набор кода, для которого он включает в себя поиск соответствующего cдля каждого a(хотя я дам значение aпозже) и zиспользование ограничения, toroot[a,c,z]а затем cобратную подстановку в окончательное выражение functionS[a,z]. У меня также есть еще одна функция, для которой есть изменение переменной functionSR[l,z]где 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];

Моя проблема в том , когда я пытаюсь найти минимум functionS[a,z]и functionSR[l,z]для некоторых aи l, скажем , a=1и l=1это дает мне ошибку. Я думаю, это связано с поведением, cкогда a=1или 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}]}

При a = 0,1 график намного более плавный.

При a = 1 на графике больше выпуклостей

Мой код плохо написан для извлечения c? Есть ли какие-то изменения, которые можно сделать? Я где-то читал, что Reduceтакже можно использовать вместо, FindRootно я все еще разбираюсь в этом. Кроме того, используется ли здесь LocalAdaptiveкак метод NIntegrate?

ОБНОВЛЕНИЕ: Обратите внимание на опечатку, я ее исправил. В сюжетах раньше я писал c=0.1и, c=1но должно быть a=0.1и a=1.

Выражения моей проблемы даются:

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

где $b=\frac{z_s}{z_h}$ и обратите внимание, что $c=c(z_s)$( c=c[z]) , Хотя в коде c=c[a,z],$c$ должен зависеть только от $z_s$( z) поскольку$a$ будет уточнено в конце.

Кроме того, возможно, есть лучший способ найти дизайн $c$. На самом деле, у меня может быть другое ограничение, где$\frac{dS}{dz_s} = 0$ (это потому, что в конце концов мне нужно свести к минимуму $S$ относительно $z_s$) и, возможно, производная от $\eqref{1}$ относительно $z_s$, чтобы их можно было использовать для поиска $c$?

Answer 1

Источник NIntegrateсообщений об ошибках можно увидеть по коэффициенту подынтегральной функции x^d/Sqrt[1-c x^d z^d], of toroot. При c > z^-3подынтегральное выражение особенное для некоторой точки области {x, 0, 1}. Более того, если NIntegrateбы можно было интегрировать через сингулярность (а с помощью, это возможно), результатом было бы комплексное число, что (предположительно) нежелательно. Чтобы продолжить, изменить переменную интегрирования к xd = x^(d+1)и применить соответствующую Methodиз здесь .

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

Кроме того, переопределите czиспользование секанса Methodи ограничьте поиск cмежду 0и z^-3.

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

(Первоначальные предположения .5 z^-3и .6 z^-3были выбраны несколько произвольно.) С этим определением czвозвращает правильное значение c, если оно существует, и z^-3вместе с FindRoot::regedсообщением об ошибке в противном случае. С этими определениями два сюжета в вопросе могут быть получены правильно следующим образом. Для a = 1,

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

Checkпредотвращает построение короткого диапазона вокруг c = 1.42, где нет решения, хотя и не устраняет соответствующие сообщения об ошибках. Второй участок, для a = .1, является

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