diff options
Diffstat (limited to 'figs/numerical.fig/simpleq/main.jl')
| -rw-r--r-- | figs/numerical.fig/simpleq/main.jl | 192 |
1 files changed, 192 insertions, 0 deletions
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() |