Initial commitv0.0

6 files changed, 430 insertions, 0 deletions

diff --git a/figs/numerical.fig/simpleq/integration.jl b/figs/numerical.fig/simpleq/integration.jl
new file mode 100644
index 0000000..0ad456a
--- /dev/null
+++ b/figs/numerical.fig/simpleq/integration.jl
@@ -0,0 +1,28 @@
+# approximate \int_a^b f using Gauss-Legendre quadratures
+function integrate_legendre(f,a,b,weights)
+  out=0
+  for i in 1:length(weights[1])
+    out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)
+  end
+  return out
+end
+# \int f*g where g is sampled at the Legendre nodes
+function integrate_legendre_sampled(f,g,a,b,weights)
+  out=0
+  for i in 1:length(weights[1])
+    out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)*g[i]
+  end
+  return out
+end
+
+
+
+# approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures
+function integrate_chebyshev(f,a,b,N)
+  out=0
+  for i in 1:N
+    out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2)
+  end
+  return out
+end
+
diff --git a/figs/numerical.fig/simpleq/iteration.jl b/figs/numerical.fig/simpleq/iteration.jl
new file mode 100644
index 0000000..f6f044a
--- /dev/null
+++ b/figs/numerical.fig/simpleq/iteration.jl
@@ -0,0 +1,55 @@
+# \hat u_n
+function hatun_iteration(e,order,d,v,maxiter)
+  # gauss legendre weights
+  weights=gausslegendre(order)
+
+  # init V and Eta
+  (V,V0,Eta,Eta0)=init_veta(weights,d,v)
+
+  # init u and rho
+  u=Array{Array{Complex{Float64}}}(undef,maxiter+1)
+  u[1]=zeros(Complex{Float64},order)
+  rho=zeros(Complex{Float64},maxiter+1)
+
+  # iterate
+  for n in 1:maxiter
+    u[n+1]=A(e,weights,Eta)\b(u[n],e,rho[n],V)
+    rho[n+1]=rhon(u[n+1],e,weights,V0,Eta0)
+  end
+
+  return (u,rho)
+end
+
+# A
+function A(e,weights,Eta)
+  N=length(weights[1])
+  out=zeros(Complex{Float64},N,N)
+  for i in 1:N
+    k=(1-weights[1][i])/(1+weights[1][i])
+    out[i,i]=k^2+4*e
+    for j in 1:N
+      y=(weights[1][j]+1)/2
+      out[i,j]+=weights[2][j]*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3)
+    end
+  end
+  return out
+end
+
+# b
+function b(u,e,rho,V)
+  out=zeros(Complex{Float64},length(V))
+  for i in 1:length(V)
+    out[i]=V[i]+2*e*rho*u[i]^2
+  end
+  return out
+end
+
+# rho_n
+function rhon(u,e,weights,V0,Eta0)
+  S=V0
+  for i in 1:length(weights[1])
+    y=(weights[1][i]+1)/2
+    S+=-weights[2][i]*(1-y)*u[i]*Eta0[i]/(2*(2*pi)^3*y^3)
+  end
+  return 2*e/S
+end
diff --git a/figs/numerical.fig/simpleq/main.jl b/figs/numerical.fig/simpleq/main.jl
new file mode 100644
index 0000000..ece0703
--- /dev/null
+++ b/figs/numerical.fig/simpleq/main.jl
@@ -0,0 +1,192 @@
+using FastGaussQuadrature
+using Printf
+using LinearAlgebra
+include("tools.jl")
+include("integration.jl")
+include("simpleq.jl")
+include("simpleq-newton.jl")
+include("iteration.jl")
+
+
+function main()
+  ## defaults
+
+  tolerance=1e-14
+  order=100
+  maxiter=21
+
+  rho=1e-6
+  e=1e-4
+
+  d=3
+  v=v_exp3d
+  a0=a0_exp3d
+
+  # plot range when plotting in x
+  xmin=0
+  xmax=100
+  nx=100
+
+  # read cli arguments
+  (params,command)=read_args(ARGS)
+
+  # read params
+  if params!=""
+    for param in split(params,";")
+      terms=split(param,"=")
+      if length(terms)!=2
+	print(stderr,"error: could not read parameter '",param,"'.\n")
+	exit(-1)
+      end
+      lhs=terms[1]
+      rhs=terms[2]
+      if lhs=="rho"
+	rho=parse(Float64,rhs)
+      elseif lhs=="e"
+	e=parse(Float64,rhs)
+      elseif lhs=="tolerance"
+	tolerance=parse(Float64,rhs)
+      elseif lhs=="order"
+	order=parse(Int64,rhs)
+      elseif lhs=="maxiter"
+	maxiter=parse(Int64,rhs)
+      elseif lhs=="xmin"
+	xmin=parse(Float64,rhs)
+      elseif lhs=="xmax"
+	xmax=parse(Float64,rhs)
+      elseif lhs=="nx"
+	nx=parse(Int64,rhs)
+      else
+	print(stderr,"unrecognized parameter '",lhs,"'.\n")
+	exit(-1)
+      end
+    end
+  end
+
+  ## run command
+  # e(rho)
+  if command=="energy_rho"
+    energy_rho(order,a0,d,v,maxiter,tolerance)
+  # d^2(rho*e(rho))
+  elseif command=="convexity"
+    ddrhoe(order,a0,d,v,maxiter,tolerance)
+  # u_n(x)
+  elseif command=="iteration"
+    iteration_ux(order,e,a0,d,v,maxiter)
+  else
+    print(stderr,"unrecognized command '",command,"'.\n")
+    exit(-1)
+  end
+
+end
+
+# read cli arguments
+function read_args(ARGS)
+  # flag
+  flag=""
+
+  # output strings
+  params=""
+  command=""
+
+  # loop over arguments
+  for arg in ARGS
+    # argument that starts with a dash
+    if arg[1]=='-'
+      # go through characters after dash
+      for char in arg[2:length(arg)]
+	
+	# set params
+	if char=='p'
+	  # raise flag
+	  flag="params"
+	else
+	  print_usage()
+	  exit(-1)
+	end
+      end
+    # arguments that do not start with a dash
+    else
+      if flag=="params"
+	params=arg
+      else
+	command=arg
+      end
+      # reset flag
+      flag=""
+    end
+  end
+
+  if command==""
+    print_usage()
+    exit(-1)
+  end
+
+  return (params,command)
+end
+
+# usage
+function print_usage()
+  print(stderr,"usage:  simpleq [-p params] <command>\n")
+end
+
+# exponential potential in 3 dimensions
+function v_exp3d(k)
+  return 8*pi/(1+k^2)^2
+end
+a0_exp3d=1.254356410591064819505409291110046864031181245635821944528
+
+# compute energy as a function of rho
+function energy_rho(order,a0,d,v,maxiter,tolerance)
+  minlrho=-6
+  maxlrho=2
+  nlrho=100
+
+  for j in 0:nlrho-1
+    rho=10^(minlrho+(maxlrho-minlrho)*j/nlrho)
+    # linear spacing
+    #rho=10.0^minlrho+(10.0^maxlrho-10.0^minlrho)*j/nlrho
+
+    (u,E)=hatu_newton(order,rho,a0,d,v,maxiter,tolerance)
+
+    @printf("% .8e % .8e\n",rho,real(E))
+
+  end
+end
+
+# compute \partial^2(e\rho) as a function of rho
+function ddrhoe(order,a0,d,v,maxiter,tolerance)
+  minlrho=-6
+  maxlrho=2
+  nlrho=100
+
+  for j in 0:nlrho-1
+    rho=10^(minlrho+(maxlrho-minlrho)*j/nlrho)
+
+    # interval
+    drho=rho*1.01
+
+    (u,E)=hatu_newton(order,rho,a0,d,v,maxiter,tolerance)
+    (up,Ep)=hatu_newton(order,rho+drho,a0,d,v,maxiter,tolerance)
+    (um,Em)=hatu_newton(order,rho-drho,a0,d,v,maxiter,tolerance)
+
+    @printf("% .8e % .8e\n",rho,real((rho+drho)*Ep+(rho-drho)*Em-2*rho*E)/drho^2)
+
+  end
+end
+
+# compute \int u_n(x) at every step
+function iteration_ux(order,e,a0,d,v,maxiter)
+  (u,rho)=hatun_iteration(e,order,d,v,maxiter)
+
+  weights=gausslegendre(order)
+
+  intun=0.
+  for n in 2:maxiter+1
+    # compute \hat u_n(0)=1/(2*rho_n)+rho_{n-1}/2*\hat u_{n-1}(0)^2
+    intun=real(1/(2*rho[n])+rho[n-1]/2*intun^2)
+    @printf("%2d % .15e\n",n-1,abs(intun-1/rho[maxiter+1]))
+  end
+end
+
+main()
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
+
diff --git a/figs/numerical.fig/simpleq/simpleq.jl b/figs/numerical.fig/simpleq/simpleq.jl
new file mode 100644
index 0000000..88209b0
--- /dev/null
+++ b/figs/numerical.fig/simpleq/simpleq.jl
@@ -0,0 +1,40 @@
+# \eta
+function eta(x,t,weights,d,v)
+  if d==2
+    return integrate_chebyshev(y->4*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y))/sqrt(((x+t)*y+abs(x-t)*(2-y))*((x+t)*(1+y)+abs(x-t)*(1-y))),0,1,length(weights))
+  elseif d==3
+    return (x>t ? 2*t/x : 2)* integrate_legendre(y->2*pi*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y)),0,1,weights)
+  end
+end
+
+# initialize V and Eta
+function init_veta(weights,d,v)
+  order=length(weights[1])
+  V=Array{Complex{Float64}}(undef,order)
+  Eta=Array{Array{Complex{Float64}}}(undef,order)
+  Eta0=Array{Complex{Float64}}(undef,order)
+  V0=v(0)
+  for i in 1:order
+    k=(1-weights[1][i])/(1+weights[1][i])
+    V[i]=v(k)
+    Eta[i]=Array{Complex{Float64}}(undef,order)
+    for j in 1:order
+      y=(weights[1][j]+1)/2
+      Eta[i][j]=eta(k,(1-y)/y,weights,d,v)
+    end
+    y=(weights[1][i]+1)/2
+    Eta0[i]=eta(0,(1-y)/y,weights,d,v)
+  end
+  return(V,V0,Eta,Eta0)
+end
+
+# inverse Fourier transform
+function u_x(x,u,weights,d)
+  order=length(weights[1])
+  if d==2
+    out=integrate_legendre_sampled(y->(1-y)/y^3*besselj(0,x*(1-y)/y)/(2*pi),u,0,1,weights)
+  elseif d==3
+    out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0,1,weights)
+  end
+  return out
+end
diff --git a/figs/numerical.fig/simpleq/tools.jl b/figs/numerical.fig/simpleq/tools.jl
new file mode 100644
index 0000000..1635a14
--- /dev/null
+++ b/figs/numerical.fig/simpleq/tools.jl
@@ -0,0 +1,20 @@
+# \Phi(x)=2*(1-sqrt(1-x))/x
+function Phi(x)
+  if abs(x)>1e-5
+    return 2*(1-sqrt(1-x))/x
+  else
+    return 1+x/4+x^2/8+5*x^3/64+7*x^4/128+21*x^5/512
+  end
+end
+# \partial\Phi
+function dPhi(x)
+  #if abs(x-1)<1e-5
+  #  @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x))
+  #end
+  if abs(x)>1e-5
+    return 1/(sqrt(1-x)*x)-2*(1-sqrt(1-x))/x^2
+  else
+    return 1/4+x/4+15*x^2/64+7*x^3/32+105*x^4/512+99*x^5/512
+  end
+end
+