1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
# \hat u(k) computed using Newton algorithm
function hatu_newton(order,rho,a0,d,v,maxiter,tolerance)
# compute gaussian quadrature weights
weights=gausslegendre(order)
# initialize V and Eta
(V,V0,Eta,Eta0)=init_veta(weights,d,v)
# initialize u, V and Eta
u=zeros(Complex{Float64},order)
for j in 1:order
# transformed k
k=(1-weights[1][j])/(1+weights[1][j])
if d==2
u[j]=2*pi*a0*rho/k
elseif d==3
u[j]=4*pi*a0*rho/k^2
end
end
# iterate
for i in 1:maxiter-1
new=u-inv(dXi(u,V,V0,Eta,Eta0,weights,rho,d))*Xi(u,V,V0,Eta,Eta0,weights,rho,d)
if(norm(new-u)/norm(u)<tolerance)
u=new
break
end
u=new
end
return (u,en(u,V0,Eta0,rho,weights,d)*rho/2)
end
# Xi(u)=0 is equivalent to the linear equation
function Xi(u,V,V0,Eta,Eta0,weights,rho,d)
order=length(weights[1])
# init
out=zeros(Complex{Float64},order)
# compute E before running the loop
E=en(u,V0,Eta0,rho,weights,d)
for i in 1:order
# k_i
k=(1-weights[1][i])/(1+weights[1][i])
# S_i
S=V[i]-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
# X_i
X=k^2/(2*E*rho)
# U_i
out[i]=u[i]-S/(2*E*(X+1))*Phi(S/(E*(X+1)^2))
end
return out
end
# derivative of Xi (for Newton)
function dXi(u,V,V0,Eta,Eta0,weights,rho,d)
order=length(weights[1])
# init
out=zeros(Complex{Float64},order,order)
# compute E before the loop
E=en(u,V0,Eta0,rho,weights,d)
for i in 1:order
# k_i
k=(1-weights[1][i])/(1+weights[1][i])
# S_i
S=V[i]-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
# X_i
X=k^2/(2*E*rho)
for j in 1:order
y=(weights[1][j]+1)/2
dS=-1/rho*(1-y)*Eta[i][j]/(2*(2*pi)^d*y^3)*weights[2][j]
dE=-1/rho*(1-y)*Eta0[j]/(2*(2*pi)^d*y^3)*weights[2][j]
dX=-k^2/(2*E^2*rho)*dE
out[i,j]=(i==j ? 1 : 0)-(dS-S*dE/E-S*dX/(X+1))/(2*E*(X+1))*Phi(S/(E*(X+1)^2))-S/(2*E^2*(X+1)^3)*(dS-S*dE/E-2*S*dX/(X+1))*dPhi(S/(E*(X+1)^2))
end
end
return out
end
# energy
function en(u,V0,Eta0,rho,weights,d)
return V0-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0,1,weights)
end
|
