yapanuwan
/
BergEC-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Revert "XZ method for new system" 

This reverts commit 6d1ffa0c0c.
This commit is contained in:
Nuwan Yapa 2024-03-08 11:54:04 -05:00
parent 8c74ac79da
commit cfd2eddad1
1 changed files with 13 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
XZ_technique.jl
View File

@ -2,14 +2,14 @@ using LinearAlgebra, Plots
include("ho_basis.jl") include("ho_basis.jl")
include("p_space.jl") include("p_space.jl")
ω = 0.5 * exp(-1im * 0.3 * pi) ω = 0.5 * exp(-1im * 0.47 * pi)
μ = 0.5 μ = 0.5
l = 0 l = 0
V1 = -8.6 V1 = -5
R1 = 1.75 R1 = sqrt(3)
V2 = 3.4 V2 = 2
R2 = 3.2 R2 = sqrt(10)
n_max = 20 n_max = 15
ns = collect(0:n_max) ns = collect(0:n_max)
ls = fill(l, n_max + 1) ls = fill(l, n_max + 1)
@ -18,34 +18,28 @@ T = sp_T_matrix(ns, ls; ω=ω, μ=μ)
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω) V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω)
n_EC = 8 n_EC = 8
train_cs = (2.2 - 0.7im) .+ (0.1 .* randn(n_EC)) .+ (0.1im * randn(n_EC)) train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC))
target_cs = range(2.0, 0.8, 6) target_cs = range(0.77, 0.22, 6)
train_E = zeros(ComplexF64, n_EC) train_E = zeros(ComplexF64, n_EC)
EC_basis = zeros(ComplexF64, (n_max + 1, length(train_cs))) EC_basis = zeros(ComplexF64, (n_max + 1, length(train_cs)))
exact_E = zeros(ComplexF64, length(target_cs)) exact_E = zeros(ComplexF64, length(target_cs))
extrapolate_E = similar(exact_E) extrapolate_E = similar(exact_E)
near_E_bound = 1.5 + 1.2im near_E = 0.2 + 0.2im
for (j, c) in enumerate(train_cs) for (j, c) in enumerate(train_cs)
H = T + c .* V H = T + c .* V
evals, evecs = eigen(H) evals, evecs = eigen(H)
i = argmin(abs.(evals .- near_E_bound)) i = argmin(abs.(evals .- near_E))
train_E[j] = evals[i] train_E[j] = evals[i]
EC_basis[:, j] = evecs[:, i] EC_basis[:, j] = evecs[:, i]
end end
N_EC = transpose(EC_basis) * EC_basis N_EC = transpose(EC_basis) * EC_basis
near_E_res = 1.5 - 0.4im
for (j, c) in enumerate(target_cs) for (j, c) in enumerate(target_cs)
mesh_p, mesh_w = get_mesh([0, 8 * exp(-0.5im)], 256) exact_E[j] = quick_pole_E((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), μ; cs_angle=0.5)
p_space_H = get_H_matrix((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), mesh_p, mesh_w, μ)
evals = eigvals(p_space_H)
i = argmin(abs.(evals .- near_E_res))
exact_E[j] = evals[i]
H = T + c .* V H = T + c .* V
H_EC = transpose(EC_basis) * H * EC_basis H_EC = transpose(EC_basis) * H * EC_basis
@ -57,5 +51,5 @@ end
scatter(real.(train_E), imag.(train_E), label="training") scatter(real.(train_E), imag.(train_E), label="training")
scatter!(real.(exact_E), imag.(exact_E), label="exact") scatter!(real.(exact_E), imag.(exact_E), label="exact")
scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated") scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
xlims!(-0.5,2.0) xlims!(-0.2,0.3)
ylims!(-1.5,1.5) ylims!(-0.3,0.3)