From e181f1cdec0d38d8a7cf1d858529320f6bd2e176 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Thu, 5 Dec 2024 18:12:55 -0500 Subject: [PATCH] Simplification --- berggren_3body_resonance.jl | 13 ++++--------- p_space_3body.jl | 13 ++++--------- p_space_3body_resonance.jl | 13 ++++--------- 3 files changed, 12 insertions(+), 27 deletions(-) diff --git a/berggren_3body_resonance.jl b/berggren_3body_resonance.jl index 37fc756..a6c46b9 100644 --- a/berggren_3body_resonance.jl +++ b/berggren_3body_resonance.jl @@ -47,18 +47,13 @@ end U = blockdiag(sparse.(U_blocks)...) end -@time "T" begin - T_blocks = [kron_sum(to_berg_basis(get_T_matrix(ks, μ), j1), to_berg_basis(get_T_matrix(ks, μ), j2)) for (j1, j2) in js] - T = blockdiag(sparse.(T_blocks)...) +@time "Block diagonal part" begin + Hb_blocks = [kron_sum(to_berg_basis(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ), j1), to_berg_basis(get_H_matrix((k, kp) -> V_l(j2, k, kp), ks, ws, μ), j2)) for (j1, j2) in js] + Hb = blockdiag(sparse.(Hb_blocks)...) end @time "T_cross" T_cross = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ws) ./ (2*μ) -@time "V1 and V2" begin - V_blocks = [kron_sum(to_berg_basis(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), j1), to_berg_basis(get_V_matrix((k, kp) -> V_l(j2, k, kp), ks, ws), j2)) for (j1, j2) in js] - V = blockdiag(sparse.(V_blocks)...) -end - E_max = 30 μω_global = 0.5 # due to simple relative coordinates @@ -72,7 +67,7 @@ E_max = 30 @time "V12_p" V12_p = W * V12_HO * transpose(W) @time "V12" V12 = transpose(U) * V12_p * U -@time "H" H = T + T_cross + V + V12 +@time "H" H = Hb + T_cross + V12 @time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1) display(evals) diff --git a/p_space_3body.jl b/p_space_3body.jl index 35a5439..1dfba2d 100644 --- a/p_space_3body.jl +++ b/p_space_3body.jl @@ -32,14 +32,9 @@ basis_size = length(js) * length(ks)^2 @assert length(basis) == basis_size "Something wrong with the basis" println("Basis size = $basis_size") -@time "T" begin - T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js] - T = blockdiag(sparse.(T_blocks)...) -end - -@time "V1" begin - V1_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js] - V1 = blockdiag(sparse.(V1_blocks)...) +@time "Block diagonal part" begin + Hb_blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1,_ ) in js] + Hb = blockdiag(sparse.(Hb_blocks)...) end E_max = 30 @@ -53,7 +48,7 @@ E_max = 30 @time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false) @time "V2" V2 = W_left * V2_HO * transpose(W_right) -@time "H" H = T + V1 + V2 +@time "H" H = Hb + V2 @time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1) display(evals) diff --git a/p_space_3body_resonance.jl b/p_space_3body_resonance.jl index c897e7a..7c9cf37 100644 --- a/p_space_3body_resonance.jl +++ b/p_space_3body_resonance.jl @@ -31,14 +31,9 @@ basis_size = length(js) * length(ks)^2 @assert length(basis) == basis_size "Something wrong with the basis" println("Basis size = $basis_size") -@time "T" begin - T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js] - T = blockdiag(sparse.(T_blocks)...) -end - -@time "V1" begin - V1_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js] - V1 = blockdiag(sparse.(V1_blocks)...) +@time "Block diagonal part" begin + Hb_blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1,_ ) in js] + Hb = blockdiag(sparse.(Hb_blocks)...) end E_max = 30 @@ -52,7 +47,7 @@ E_max = 30 @time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false) @time "V2" V2 = W_left * V2_HO * transpose(W_right) -@time "H" H = T + V1 + V2 +@time "H" H = Hb + V2 @time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1) display(evals)