Skip to main content

numerical integration - Elliptic integrals


In trying to reproduce results from one paper I stumbled upon a problem with definition of some elliptic integrals (this is my guess of what could be the problem).


I will first present in a simplified form what I am trying to calculate, details are in the original paper (PRL 99, 226801, see Google Scholar for PDF)





The goal is to compute the following 2D integral


I(k)=−∬Ωcd2→q4πq(1−cosθ(→k,→k−→q)),


where θ(→a,→b) is the angle between →a and →b, q=|→q|. For those who would like to compare with the paper, this is essentially Eq.(2a), where for simplicity I set e=κ=1, selected the case s=1 and substituted all the definitions into one equation.


The integration domain is Ωc:|k|≤kc, where kc is a positive number.


The analytic result is known to be [cf. Eq. (3a)]:


I(k)=1πkc[h(k/kc)−f(k/kc)],I(0)=−12kc.


Assuming, we want to know the result for $k, there are following definitions


f(x)=E(x),h(x)=x[π4log(4/x)−π8]−x∫x0dyy−3[K(y)−E(y)−π4y2].

Here K(x) and E(x) are the complete elliptic integral of the first and second kinds, respectively. I do not know, how this integral can be computed, neither by hand or with mathematica...


What is disturbing is that I was not able verify the integral numerically.


In the following, I will first rewrite all the equations in MA language.



i[1]=Integrate[EllipticK[y^2]-EllipticE[y^2],{y,0,1/x},
Assumptions->x>1]
i[2]=Integrate[(EllipticK[y^2]-EllipticE[y^2]-Ï€/4 y^2)/y^3,{y,0,x},
Assumptions->0<=x<=1]
f[x_]=EllipticE[x^2]
h[x_]=x(Ï€/4Log[4/x]-Ï€/8)-x i[2]

Notice, it takes a while to compute i[2]. Now, we will be interested in the ΔI(k)=I(k)−I(0) function


xi[k_,kc_]:=kc/Ï€(h[k/kc]-f[k/kc])
Δxi[k_,kc_]:=xi[k,kc]+kc/2


Now we define the numerical integral (adding a small cut-off a) transforming it to polar coordinates and assuming →k∥→ex


Δni[k_?NumericQ,kc_?NumericQ,a_?NumericQ]:=1/(4 π) NIntegrate[((k- q Cos[θ])/Sqrt[k^2+q^2-2 k q Cos[θ]]),{q,a,kc},{θ,0,2π},PrecisionGoal->4]

and compare


dataI=Table[{k,Δni[k,30,10^-5]},{k,0.1,2,0.1}]
Plot[Δxi[k,30],{k,0,2},Epilog->{PointSize[Medium],Point[dataI]},PlotRange->{0,2.3}]

enter image description here


The points should exactly fall onto the analytic curve, but they are not... I would be happy with any of the answers:




  1. showing how the integral can be computed analytically with MA starting from the definition (basically confirming that my interpretation of the formula in the paper is correct),

  2. fixing a problem with MA numerics.


Notice, I can easily verify Fig.1 of that paper with MA. But the integral considered here is not plotted there.



Answer



The article "Density Dependent Exchange Contribution to ∂μ/∂n and Compressibility in Graphene" by E.H. Hwang,Ben Yu-Kuang Hu, and S. Das Sarma has a typo in the definition of h (there must be a plus before π8). After correction, the results match (I wrote down the finished results for the integrals, so as not to waste time every time to calculate them)


(*i[1]=Integrate[EllipticK[y^2]-EllipticE[y^2],{y,0,1/x},Assumptions\
\[Rule]x>1]
i[2]=Integrate[(EllipticK[y^2]-EllipticE[y^2]-\[Pi]/4 \

y^2)/y^3,{y,0,x},Assumptions\[Rule]0\[LessEqual]x\[LessEqual]1]*)

i1[x_] := (\[Pi] (-HypergeometricPFQ[{-(1/2), 1/2, 1/2}, {1, 3/2}, 1/
x^2] + HypergeometricPFQ[{1/2, 1/2, 1/2}, {1, 3/2}, 1/x^2]))/(
2 x)
i2[x_] :=
3/256 \[Pi] x^2 (HypergeometricPFQ[{1, 1, 3/2, 5/2}, {2, 3, 3},
x^2] + 3 HypergeometricPFQ[{1, 1, 5/2, 5/2}, {2, 3, 3}, x^2])
f[x_] := If[x <= 1, EllipticE[x^2],
x EllipticE[1/x^2] - (x - 1/x) EllipticK[1/x^2]]

h[x_] := If[x <= 1, x (\[Pi]/4 Log[4/x] + \[Pi]/8) - x i2[x], x i1[x]]

xi[k_, kc_] := kc/\[Pi] (h[k/kc] - f[k/kc])
\[CapitalDelta]xi[k_, kc_] := xi[k, kc] + kc/2

\[CapitalDelta]ni[k_?NumericQ, kc_?NumericQ, a_?NumericQ] :=
1/(4 \[Pi]) NIntegrate[((k - q Cos[\[Theta]])/
Sqrt[k^2 + q^2 - 2 k q Cos[\[Theta]]]), {q, a, kc}, {\[Theta], 0,
2 \[Pi]}]


dataI = Table[{k, \[CapitalDelta]ni[k, 30, 10^-10]}, {k, 0.1, 2, 0.1}]
Plot[\[CapitalDelta]xi[k, 30], {k, 0, 2},
Epilog -> {PointSize[Medium], Point[dataI]}, PlotRange -> {0, 2.3}]

Figure 1


Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]