diff options
Diffstat (limited to 'figs/numerical.fig')
| -rw-r--r-- | figs/numerical.fig/Makefile | 32 | ||||
| -rw-r--r-- | figs/numerical.fig/bosegas.gnuplot | 34 | ||||
| -rw-r--r-- | figs/numerical.fig/convexity.gnuplot | 30 | ||||
| -rw-r--r-- | figs/numerical.fig/holzmann_2019-09-28.dat | 9 | ||||
| -rw-r--r-- | figs/numerical.fig/simpleq/integration.jl | 28 | ||||
| -rw-r--r-- | figs/numerical.fig/simpleq/iteration.jl | 55 | ||||
| -rw-r--r-- | figs/numerical.fig/simpleq/main.jl | 192 | ||||
| -rw-r--r-- | figs/numerical.fig/simpleq/simpleq-newton.jl | 95 | ||||
| -rw-r--r-- | figs/numerical.fig/simpleq/simpleq.jl | 40 | ||||
| -rw-r--r-- | figs/numerical.fig/simpleq/tools.jl | 20 |
10 files changed, 535 insertions, 0 deletions
diff --git a/figs/numerical.fig/Makefile b/figs/numerical.fig/Makefile new file mode 100644 index 0000000..8fb4411 --- /dev/null +++ b/figs/numerical.fig/Makefile @@ -0,0 +1,32 @@ +PROJECTNAME=bosegas convexity + +DATS=erho.dat ddrhoe.dat +PDFS=$(addsuffix .pdf, $(PROJECTNAME)) +TEXS=$(addsuffix .tikz.tex, $(PROJECTNAME)) + +all: $(PDFS) + +$(PDFS): $(DATS) + gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@) + pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@) + +erho.dat: + julia ./simpleq/main.jl -p "tolerance=1e-14;order=100" energy_rho > $@ +ddrhoe.dat: + julia ./simpleq/main.jl -p "tolerance=1e-14;order=100" convexity > $@ + +install: $(PDFS) + cp $^ $(INSTALLDIR)/ + +clean-aux: + rm -f $(addsuffix .tikz.tex, $(PROJECTNAME)) + rm -f $(addsuffix .aux, $(PROJECTNAME)) + rm -f $(addsuffix .log, $(PROJECTNAME)) + +clean-dat: + rm -f $(DATS) + +clean-tex: + rm -f $(PDFS) + +clean: clean-aux clean-tex diff --git a/figs/numerical.fig/bosegas.gnuplot b/figs/numerical.fig/bosegas.gnuplot new file mode 100644 index 0000000..bfcffac --- /dev/null +++ b/figs/numerical.fig/bosegas.gnuplot @@ -0,0 +1,34 @@ +set ylabel "$\\displaystyle\\frac{e}{4\\pi\\rho}$" norotate offset -1,0 +set xlabel "$\\rho$" + +set xtics 1e-6, 100, 100 +set xtics add ("$10^{-6}$" 0.000001, "$10^{-4}$" 0.0001, "$10^{-2}$" 0.01, "$1$" 1.0, "$10^2$" 100) +unset mxtics + +set ytics 0.6, 0.1 +set mytics 2 + +set xrange [0.000001:100] +set yrange [0.6:1.05] + +# default output canvas size: 12.5cm x 8.75cm +set term lua tikz size 8,6 standalone + +set key off + + +# set linestyle +set style line 1 linetype rgbcolor "#4169E1" linewidth 3 +set style line 2 linetype rgbcolor "#DC143C" linewidth 3 +set style line 3 linetype rgbcolor "#32CD32" linewidth 3 +set style line 4 linetype rgbcolor "#4B0082" linewidth 3 +set style line 5 linetype rgbcolor "#DAA520" linewidth 3 + +set pointsize 1 + +set logscale x + +plot "erho.dat" using 1:(0.95*$2/$1/(4*pi)):(1.05*$2/$1/(4*pi)) with filledcurves linetype rgbcolor "#DDDDDD", \ + "erho.dat" using 1:($2/$1/(4*pi)) with lines linestyle 1, \ + "holzmann_2019-09-28.dat" using 1:($2/$1/(4*pi)) with points linestyle 2 + diff --git a/figs/numerical.fig/convexity.gnuplot b/figs/numerical.fig/convexity.gnuplot new file mode 100644 index 0000000..f6f606b --- /dev/null +++ b/figs/numerical.fig/convexity.gnuplot @@ -0,0 +1,30 @@ +set ylabel "$\\frac1{4\\pi}\\partial_\\rho^2(e\\rho)$" offset -4,0 norotate +set xlabel "$\\rho$" + +set xtics 1e-6, 100, 100 +set xtics add ("$10^{-6}$" 0.000001, "$10^{-4}$" 0.0001, "$10^{-2}$" 0.01, "$1$" 1.0, "$10^2$" 100) +unset mxtics +set xrange [0.000001:100] + +set ytics 1.2, 0.2 +set mytics 2 +set yrange [1.2:2.1] + +# default output canvas size: 12.5cm x 8.75cm +set term lua tikz size 8,6 standalone + +# no key +set key off + +# set linestyle +set style line 1 linetype rgbcolor "#4169E1" linewidth 3 +set style line 2 linetype rgbcolor "#DC143C" linewidth 3 +set style line 3 linetype rgbcolor "#000000" linewidth 1 +set style line 4 linetype rgbcolor "#4B0082" linewidth 3 +set style line 5 linetype rgbcolor "#DAA520" linewidth 3 + +set pointsize 0.6 + +set logscale x + +plot "ddrhoe.dat" using 1:($2/(4*pi)) with lines linestyle 1 diff --git a/figs/numerical.fig/holzmann_2019-09-28.dat b/figs/numerical.fig/holzmann_2019-09-28.dat new file mode 100644 index 0000000..fe9c0e5 --- /dev/null +++ b/figs/numerical.fig/holzmann_2019-09-28.dat @@ -0,0 +1,9 @@ +## data from M. Holzmann, 2019-09-22 +# rho E0 E0+dE0 +1e-4 8.3500e-4 8.3600e-4 +1e-3 9.1340e-3 9.1350e-3 +1e-2 1.0609e-1 1.0610e-1 +1e-1 1.1930e+0 1.1940e+0 +1e-0 1.2445e+1 1.2446e+1 +1e+1 1.2544e+2 1.2545e+2 + 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 + |