1
$\begingroup$

I want to solve the following problem.

First, the setup. A function $V(\varphi)$ defined as \begin{multline} V(\varphi)=\frac{\alpha-3}{4}e^{\sqrt{2\alpha}\varphi}+\frac{\alpha-5} {\gamma}e^{(\alpha-1)\sqrt{\frac{2}{\alpha}}\varphi}+\frac{\alpha^2- 7\alpha+4}{\alpha\gamma^2}e^{(\alpha-2)\sqrt{\frac{2}{\alpha}}\varphi}-\\ -\frac{\alpha(\gamma+2)^2}{4\gamma^2}+\frac{(\gamma+2)(3\gamma+14)} {4\gamma^2}-\frac{4}{\alpha\gamma^2}~, \end{multline} where $\varphi$ is a real variable and $\alpha$ is a constant real parameter (specifically I'm interested in $\alpha>(7+\sqrt{33})/2$). $\gamma$ is just $\gamma=2(2-\alpha)/\alpha$.

Then, there is the integral: $$N_e=\int^{\varphi_i}_{\varphi_f} d \varphi\frac{V(\varphi)}{V'(\varphi)},$$ where $\varphi_i$ is defined as $$\varphi_i=\sqrt{\frac{\alpha}{2}}\ln\left(\frac{\alpha^2-7\alpha+4}{(\alpha-3)(\alpha-2)}\right)~,$$ and is negative for $\alpha>(7+\sqrt{33})/2$, while $\varphi_f$ is the (negative) root of the equation $$\frac{1}{2}\left(\frac{V'(\varphi)}{V(\varphi)}\right)^2=1~.$$

My goal is to find the value of $\alpha$ for which $N_e=50$.

My attempt:

Clear["Global`*"]

V[φ_, α_] := ((α - 3)/4)*E^(Sqrt[2*α]*φ) + 
    ((α - 5)/γ)*E^((α - 1)*Sqrt[2/α]*φ) + ((α^2 - 7*α + 4)/(α*γ^2))*
     E^((α - 2)*Sqrt[2/α]*φ) - (α*(γ + 2)^2)/(4*γ^2) + 
    ((γ + 2)*(3*γ + 14))/(4*γ^2) - 4/(α*γ^2); 

γ = (-2*(α - 2))/α; 

ϵ[φ_, α_] := (1/2)*(D[V[φ, α], φ]/V[φ, α])^2; 

φi[α_] := Sqrt[α/2]*Log[(α^2 - 7*α + 4)/((α - 3)*(α - 2))]; 

φf[α_] := φ /. FindRoot[ϵ[φ, α] == 1, {φ, -0.3}]; 

Ne[(α_)?NumericQ] := NIntegrate[V[φ, α]/D[V[φ, α], φ], {φ, φf[α], φi[α]}];

FindRoot[Ne[α] == 50, {α, 11}]

However, this gives a list of errors most of which are on NIntegrate. For example the first two are

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in φ near {φ} = {-0.748946837826036636752475288418998372825970122595301559087488385558}. NIntegrate obtained 25.416144461679007` and 2.053555362055963` for the integral and error estimates.

And the last one on FindRoot:

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

I'm not sure how to improve my code. Desirable accuracy for $\alpha$ would be two digits after the decimal point.

Updated the code.

$\endgroup$
6
  • $\begingroup$ Please check the definition of \[Epsilon][\[Alpha]_]=... (rigth hand side only depends on \[CurlyPhi] $\endgroup$ Jun 12, 2019 at 10:16
  • $\begingroup$ @UlrichNeumann Should I define it as [Epsilon][[Alpha]_,[CurlyPhi]] ? Because V contains $\alpha$ as a parameter. Or should I define V[[CurlyPhi]_,[Alpha]] too? $\endgroup$
    – Kosm
    Jun 12, 2019 at 11:03
  • $\begingroup$ @ Kosm I think you should define both functions depending on \[CurlyPhi]_, \[Alpha]_ ! $\endgroup$ Jun 12, 2019 at 11:06
  • $\begingroup$ @UlrichNeumann Updated my code. $\endgroup$
    – Kosm
    Jun 12, 2019 at 11:55
  • $\begingroup$ @ Kosm: Reflecting about your problem, are you sure about the integrand V/V'. What about V'/v? $\endgroup$ Jun 12, 2019 at 18:46

2 Answers 2

4
+50
$\begingroup$

Your integrand has a simple pole at φi[α] and therefore does not converge.

First fix up your definitions. If you want γ defined outside of V[], define it with complete arguments γ[α_] :=....

V[φ_, α_] := With[{γ = (-2*(α - 2))/α},
   ((α - 3)/4)*E^(Sqrt[2*α]*φ)
     + ((α - 5)/γ)*E^((α - 1)*Sqrt[2/α]*φ)
     + ((α^2 - 7*α + 4)/(α*γ^2))*E^((α - 2)*Sqrt[2/α]*φ)
     - (α*(γ + 2)^2)/(4*γ^2)
     + ((γ + 2)*(3*γ + 14))/(4*γ^2)
     - 4/(α*γ^2)
   ];

ϵ[φ_, α_] := (1/2)*(D[V[φ, α], φ]/V[φ, α])^2;

Here we can see that the integrand has a pole if the coefficient does not vanish (which it does not seem to do):

Assuming[α > 0,
 Simplify@ SeriesCoefficient[V[φ, α]/D[V[φ, α], φ], {φ, φi[α], -1}]
 ]
(*
(12 - α^3 - 
 12 ((4 - 7 α + α^2)/(6 - 5 α + α^2))^-α + α^2 (7 - 
    3 ((4 - 7 α + α^2)/(6 - 5 α + α^2))^-α) + α (-16 + 
    21 ((4 - 7 α + α^2)/(6 - 5 α + α^2))^-α))/((-2 + α)^2 (-3 - 2 α + α^2))
*)
$\endgroup$
4
  • $\begingroup$ I suppose that the pole with the largest value of $\alpha$ is at $\alpha=(7+\sqrt{33})/2$. Since I'm interested in $\alpha$ greater than this value, should I add the assumption $\alpha>(7+\sqrt{33})/2$ to NIntegrate? $\endgroup$
    – Kosm
    Jun 13, 2019 at 19:43
  • 1
    $\begingroup$ @Kosm You have the pole at φ == φi[α] as a limit of integration in NIntegrate. I don't see how the assumption helps with the computation. Besides, NIntegrate never sees the variable α because a numeric value is always substituted for. (If the interval of integration contained both sides of the pole, you could use the principal value; if it did not contain the pole, the integral would converge. Aside from there being some such error in the setup of the integral, there's nothing I can suggest at present.) $\endgroup$
    – Michael E2
    Jun 13, 2019 at 20:44
  • $\begingroup$ @Kosm I see you've put a bounty on the question. The answer "the integral does not converge" (that is, it diverges) is a complete mathematical answer. One might say it this way in this case, "the integral goes to infinity for all α." This lack of mathematical convergence is a distinct concept from the lack of numerical convergence reported by NIntegrate[], although mathematical divergence implies numerical divergence. In other words, if the answer cannot be infinity for other reasons, then the integral has been formulated incorrectly. $\endgroup$
    – Michael E2
    Jun 18, 2019 at 15:23
  • $\begingroup$ I see. I guess I can shift the $\varphi_i$ by a small amount away from the pole then. $\endgroup$
    – Kosm
    Jun 19, 2019 at 8:06
1
$\begingroup$

Here my solution approach:

Clear["Global`*"]
\[Gamma] = (-2*(\[Alpha] - 2))/\[Alpha];
V[\[CurlyPhi]_, \[Alpha]_] := ((\[Alpha] - 3)/4)*E^(Sqrt[2*\[Alpha]]*\[CurlyPhi]) + ((\[Alpha] - 5)/\[Gamma])*E^((\[Alpha] - 1)*Sqrt[2/\[Alpha]]*\[CurlyPhi]) + ((\[Alpha]^2 - 7*\[Alpha] +4)/(\[Alpha]*\[Gamma]^2))*E^((\[Alpha] - 2)*Sqrt[2/\[Alpha]]*\[CurlyPhi]) - (\[Alpha]*(\[Gamma] +2)^2)/(4*\[Gamma]^2) + ((\[Gamma] + 2)*(3*\[Gamma] +14))/(4*\[Gamma]^2) - 4/(\[Alpha]*\[Gamma]^2);
\[CurlyPhi]i[\[Alpha]_] :=Sqrt[\[Alpha]/2]*Log[(\[Alpha]^2 - 7*\[Alpha] + 4)/((\[Alpha] - 3)*(\[Alpha] - 2))];
\[Epsilon][\[CurlyPhi]_?NumericQ, \[Alpha]_?NumericQ] := (1/2)(Derivative[1, 0][V][\[CurlyPhi], \[Alpha]]/V[\[CurlyPhi], \[Alpha]])^2;

The solution \[CurlyPhi]f[\[Alpha]] doesn't work in the final evaluation, taht 's why I try to find a good approximation

cplot = ContourPlot[\[Epsilon][\[CurlyPhi], \[Alpha]] ==1, {\[Alpha], (7 + Sqrt[33])/2., (7 + Sqrt[33])/2. +Pi}, {\[CurlyPhi], -Pi/2, -.010}, FrameLabel -> {\[Alpha], \[CurlyPhi]}];
\[Alpha]\[CurlyPhi] = cplot[[1, 1]]; (*plotpoints*)
\[CurlyPhi]f[\[Alpha]_] := Fit[\[Alpha]\[CurlyPhi], {1, \[Alpha], \[Alpha]^2},\[Alpha]] //Evaluate
Show[{cplot,Plot[\[CurlyPhi]f[\[Alpha]], {\[Alpha], 6, 10},PlotStyle ->{Opacity[.1], Thickness[.02], Red}]}, PlotRange -> {-1, 0}]

enter image description here

Now it's possible to define function Ne[\[Alpha]].

Ne[(\[Alpha]_)?NumericQ] := Evaluate@Apply[NIntegrate, {V[\[CurlyPhi], \[Alpha]]/Derivative[1, 0][V][\[CurlyPhi],\[Alpha]], {\[CurlyPhi],\[CurlyPhi]f[\[Alpha]], \[CurlyPhi]i[\[Alpha]]}}];

Plot[Ne[\[Alpha]], {\[Alpha], 6.5, 15}, Evaluated -> True,PlotRange -> {0, 100} ,GridLines -> {{11}, {50}}, MaxRecursion -> 3,AxesLabel -> {\[Alpha], Ne[\[Alpha]]}]

enter image description here

FindRoot doesn't solve the equation Ne[\[Alpha]] == 50 but

NMinimize[{1,Ne[\[Alpha]] == 50, \[Alpha] > (7 + Sqrt[33])/2}, \[Alpha]]
(*{1., {\[Alpha] -> 10.8255}}*)

does!

$\endgroup$
10
  • $\begingroup$ This is strange. It gives me again a list of errors starting with "NIntegrate::nlim: [CurlyPhi] = -1.10127+0.0972841 [Alpha]-0.00359968 [Alpha]^2 is not a valid limit of integration." $\endgroup$
    – Kosm
    Jun 13, 2019 at 5:31
  • $\begingroup$ NMinimize gives this: NMinimize::nosat: Obtained solution does not satisfy the following constraints within Tolerance -> 0.001`: {50-Ne[[Alpha]]==0}. $\endgroup$
    – Kosm
    Jun 13, 2019 at 5:31
  • $\begingroup$ It still solves NMinimize, but on my machine it is "[Alpha] -> 11.17" $\endgroup$
    – Kosm
    Jun 13, 2019 at 5:32
  • $\begingroup$ I edited my answer Ne[] now without message "valid limit of integration". Convergence problem of NIntegrate still remains! $\endgroup$ Jun 13, 2019 at 6:32
  • $\begingroup$ why do you think the solutions differ for us? $\endgroup$
    – Kosm
    Jun 13, 2019 at 6:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.