From e6bf8349d7fc554bdbd915cc65c44f23c1b86e75 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian@jauslin.org>
Date: Mon, 16 Dec 2019 15:06:59 -0500
Subject: Initial commit

---
 figs/numerical.fig/simpleq/main.jl | 192 +++++++++++++++++++++++++++++++++++++
 1 file changed, 192 insertions(+)
 create mode 100644 figs/numerical.fig/simpleq/main.jl

(limited to 'figs/numerical.fig/simpleq/main.jl')

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()
-- 
cgit v1.2.3-70-g09d2