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

Notebooks -> Scripts

This commit is contained in:
Nuwan Yapa 2024-12-19 17:07:19 -05:00
parent b3934d6885
commit 9c83e3e56b
6 changed files with 99 additions and 232 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

100
EC_test.ipynb
View File

@ -1,100 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"include(\"p_space.jl\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"vertices = [0, 0.4 - 0.15im, 0.8, 6]\n",
"subdivisions = [128, 128, 128]\n",
"p, w = get_mesh(vertices, subdivisions)\n",
"mesh_E = p.*p ./ (2*0.5)\n",
"\n",
"# ResonanceEC: Eq. (20)\n",
"V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"training_points = range(0.75, 0.45, 5)\n",
"training_E = Vector{ComplexF64}(undef, length(training_points))\n",
"EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))\n",
"\n",
"for (j, c) in enumerate(training_points)\n",
" evals, evecs = eigen(get_H_matrix(V_system(c), p, w))\n",
" i = identify_pole_i(p, evals)\n",
" training_E[j] = evals[i]\n",
" EC_basis[:, j] = evecs[:, i]\n",
"end\n",
"\n",
"scatter(real.(training_E), imag.(training_E), label=\"training\")\n",
"plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n",
"xlims!(0,1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"extrapolate_points = range(0.40, 0.25, 5)\n",
"ref_E = 0.2 - 0.1im\n",
"\n",
"exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n",
"extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n",
"\n",
"EC_basis_w = EC_basis .* w\n",
"N_EC = transpose(EC_basis_w) * EC_basis\n",
"\n",
"for (j, c) in enumerate(extrapolate_points)\n",
" exact_E[j] = quick_pole_E(V_system(c))\n",
"\n",
" EC_basis_w = EC_basis .* w\n",
" H = get_H_matrix(V_system(c), p, w)\n",
" H_EC = transpose(EC_basis_w) * H * EC_basis\n",
" evals = eigvals(H_EC, N_EC)\n",
" i = argmin(abs.(evals .- ref_E))\n",
" ref_E = evals[i]\n",
" extrapolate_E[j] = evals[i]\n",
"end\n",
"\n",
"scatter(real.(training_E), imag.(training_E), label=\"training\")\n",
"scatter!(real.(exact_E), imag.(exact_E), label=\"exact\")\n",
"scatter!(real.(extrapolate_E), imag.(extrapolate_E), label=\"extrapolated\")\n",
"plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n",
"xlims!(0,1)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.10.2",
"language": "julia",
"name": "julia-1.10"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

48
calculations/R2R_Berggren.jl Normal file
View File

@ -0,0 +1,48 @@
using Plots
include("../p_space.jl")
# contour
vertices = [0, 0.4 - 0.15im, 0.8, 6]
subdivisions = [128, 128, 128]
p, w = get_mesh(vertices, subdivisions)
mesh_E = p.*p ./ (2*0.5)
# ResonanceEC: Eq. (20)
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
training_points = range(0.75, 0.45, 5)
training_E = Vector{ComplexF64}(undef, length(training_points))
EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
for (j, c) in enumerate(training_points)
evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
i = identify_pole_i(p, evals)
training_E[j] = evals[i]
EC_basis[:, j] = evecs[:, i]
end
extrapolate_points = range(0.40, 0.25, 5)
ref_E = 0.2 - 0.1im
exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
EC_basis_w = EC_basis .* w
N_EC = transpose(EC_basis_w) * EC_basis
for (j, c) in enumerate(extrapolate_points)
exact_E[j] = quick_pole_E(V_system(c))
H = get_H_matrix(V_system(c), p, w)
H_EC = transpose(EC_basis_w) * H * EC_basis
evals = eigvals(H_EC, N_EC)
i = argmin(abs.(evals .- ref_E))
global ref_E = evals[i]
extrapolate_E[j] = evals[i]
end
scatter(real.(training_E), imag.(training_E), label="training")
scatter!(real.(exact_E), imag.(exact_E), label="exact")
scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
plot!(real.(mesh_E), imag.(mesh_E), label="contour")
xlims!(0,1)

62
p_space_test.ipynb
View File

@ -1,62 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"include(\"helper.jl\")\n",
"include(\"p_space.jl\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"vertices = [0, 1 - 0.5im, 2, 6]\n",
"subdivisions = [64, 64, 64]\n",
"p, w = get_mesh(vertices, subdivisions)\n",
"\n",
"scatter(real.(p), imag.(p))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ComplexScaling-FV: Eq. (54)\n",
"V_pq(p, q) = -10 * g1(1, p, q)\n",
"\n",
"H = get_H_matrix(V_pq, p, w)\n",
"evals = eigen(H).values\n",
"\n",
"E_target = 0.258 - 0.164im\n",
"E = nearest(evals, E_target)\n",
"\n",
"print(\"E = $E\")\n",
"scatter(real.(evals), imag.(evals), xlim = (0,1))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.9.0",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

70
simple_berggren.ipynb
View File

@ -1,70 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"include(\"helper.jl\")\n",
"include(\"p_space.jl\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"vertices = [0, 0.5 - 0.3im, 1, 6]\n",
"subdivisions = [64, 64, 64]\n",
"p, w = get_mesh(vertices, subdivisions)\n",
"\n",
"# resonance example from my thesis\n",
"V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)\n",
"\n",
"basis_eig = eigen(get_H_matrix(V_basis, p, w))\n",
"basis = basis_eig.vectors .* w\n",
"\n",
"evals = basis_eig.values\n",
"E_target = 0.7\n",
"E = include(\"helper.jl\")\n",
"\n",
"print(\"E = $E\")\n",
"scatter(real.(evals), imag.(evals), xlim = (0,2))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ResonanceEC: Eq. (20)\n",
"V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))\n",
"\n",
"H = get_H_matrix(V_system(0.45), p, w)\n",
"H_berggren = transpose(basis) * H * basis\n",
"\n",
"evals = eigvals(H)\n",
"scatter(real.(evals), imag.(evals), xlim = (0, 0.5))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.10.2",
"language": "julia",
"name": "julia-1.10"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

21
test/p_space_basic.jl Normal file
View File

@ -0,0 +1,21 @@
using Plots
include("../helper.jl")
include("../p_space.jl")
vertices = [0, 1 - 0.5im, 2, 6]
subdivisions = [64, 64, 64]
p, w = get_mesh(vertices, subdivisions)
scatter(real.(p), imag.(p))
# ComplexScaling-FV: Eq. (54)
V_pq(p, q) = -10 * g1(1, p, q)
H = get_H_matrix(V_pq, p, w)
evals = eigen(H).values
E_target = 0.258 - 0.164im
E = nearest(evals, E_target)
print("E = $E")
scatter(real.(evals), imag.(evals), xlim = (0,1))

30
test/simple_berggren.jl Normal file
View File

@ -0,0 +1,30 @@
using Plots
include("../helper.jl")
include("../p_space.jl")
vertices = [0, 0.5 - 0.3im, 1, 6]
subdivisions = [64, 64, 64]
p, w = get_mesh(vertices, subdivisions)
# resonance example from my thesis
V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)
basis_eig = eigen(get_H_matrix(V_basis, p, w))
basis = basis_eig.vectors .* w
basis_evals = basis_eig.values
E_target = 0.7
E = nearest(basis_evals, E_target)
print("pole E = $E")
# ResonanceEC: Eq. (20)
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
H = get_H_matrix(V_system(0.45), p, w)
H_berggren = transpose(basis) * H * basis
evals = eigvals(H)
scatter(real.(evals), imag.(evals), label="E")
scatter!(real.(basis_evals), imag.(basis_evals), label="basis", markershape=:x)
xlims!((0, 1))