Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat (limited to 'figs/numerical.fig/simpleq/simpleq-newton.jl')
-rw-r--r--figs/numerical.fig/simpleq/simpleq-newton.jl95
1 files changed, 95 insertions, 0 deletions
diff --git a/figs/numerical.fig/simpleq/simpleq-newton.jl b/figs/numerical.fig/simpleq/simpleq-newton.jl
new file mode 100644
index 0000000..2c70162
--- /dev/null
+++ b/figs/numerical.fig/simpleq/simpleq-newton.jl
@@ -0,0 +1,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
+