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/simpleq-newton.jl | 95 ++++++++++++++++++++++++++++
 1 file changed, 95 insertions(+)
 create mode 100644 figs/numerical.fig/simpleq/simpleq-newton.jl

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

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