Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Callback finds only one solution #1040

Open
IromainI opened this issue Jul 15, 2024 · 1 comment
Open

Callback finds only one solution #1040

IromainI opened this issue Jul 15, 2024 · 1 comment
Labels

Comments

@IromainI
Copy link

Hi everyone!
I am trying to find the right parameters that satisfy the solved ODE system. In my problem, it is necessary that the solution r[5]. For verification, I use I use callback return abs(real(u[5])) >= 0.1. When this condition is met, I want the current solution to be interrupted and a new one with new parameters to start. The problem is as follows: when the first solution is found, the program turns off, the function affect! is not called. Thank you

These are the initial values of the variables

using DifferentialEquations

global_index_t1 = 1
global_index_t2 = 1
global_index_Ω1 = 1
global_index_Ω2 = 1
global_index_Δ1 = 1
global_index_Δ2 = 1
global_index_Γ1 = 1
global_index_Γ2 = 1
global_index_γ1 = 1
global_index_γ2 = 1
global_index_γ3 = 1
global_index_σ1 = 1
global_index_σ2 = 1
finish = false

successful_combinations = []

Ω = [10.0, 10.0]
Δ = [0.0, 0.0]
Γ = [0.5, 0.5]
γ = [1/4, 0.5, 1/4]
t0 = [24.0, 20.0]
σ = [3.0, 3.0]

Ω01 = collect(9:1:11)
Ω02 = collect(9:1:11)
Δ01 = collect(0.0:1:2.0)
Δ02 = collect(0.0:1:1.0)
σ01 = collect(2.0:1:4.0)
σ02 = collect(2.0:1:4.0)
t00 = collect(23:1:25)
t01 = collect(18:1:22)
Γ01 = collect(0.4:0.1:0.6)
Γ02 = collect(0.4:0.1:0.6)
γ01 = collect(0.2:0.05:0.3)
γ02 = collect(0.0:0.05:0.1)
γ03 = collect(0.2:0.05:0.3)

This is the part with callack

function condition(u, t, integrator)
    return abs(real(u[5])) >= 0.1
end

function affect!(integrator)
    global global_index_t1, global_index_t2, global_index_Ω1, global_index_Ω2, global_index_Δ1, global_index_Δ2
    global global_index_Γ1, global_index_Γ2, global_index_γ1, global_index_γ2, global_index_γ3
    global global_index_σ1, global_index_σ2, finish

    println(#"Event triggered at t = ", integrator.t,
            global_index_t1,
            global_index_t2,
            global_index_Ω1,
            global_index_Ω2,
            global_index_Δ1,
            global_index_Δ2,
            global_index_Γ1,
            global_index_Γ2,
            global_index_γ1,
            global_index_γ2,
            global_index_γ3,
            global_index_σ1,
            global_index_σ2)

            if global_index_t1 < length(t00)
                global_index_t1 += 1
            elseif global_index_t2 < length(t01)
                global_index_t1 = 1
                global_index_t2 += 1
            elseif global_index_Ω1 < length(Ω01)
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_Ω1 += 1
            elseif global_index_Ω2 < length(Ω02)
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_Ω2 += 1
            elseif global_index_Δ1 < length(Δ01)
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_Δ1 += 1
            elseif global_index_Δ2 < length(Δ02)
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_Δ2 += 1
            elseif global_index_Γ1 < length(Γ01)
                global_index_Δ2 = 1
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_Γ1 += 1
            elseif global_index_Γ2 < length(Γ02)
                global_index_Γ1 = 1
                global_index_Δ2 = 1
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_Γ2 += 1
            elseif global_index_γ1 < length(γ01)
                global_index_Γ2 = 1
                global_index_Γ1 = 1
                global_index_Δ2 = 1
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_γ1 += 1
            elseif global_index_γ2 < length(γ02)
                global_index_γ1 = 1
                global_index_Γ2 = 1
                global_index_Γ1 = 1
                global_index_Δ2 = 1
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_γ2 += 1
            elseif global_index_γ3 < length(γ03)
                global_index_γ2 = 1
                global_index_γ1 = 1
                global_index_Γ2 = 1
                global_index_Γ1 = 1
                global_index_Δ2 = 1
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_γ3 += 1
            elseif global_index_σ1 < length(σ01)
                global_index_γ3 = 1
                global_index_γ2 = 1
                global_index_γ1 = 1
                global_index_Γ2 = 1
                global_index_Γ1 = 1
                global_index_Δ2 = 1
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_σ1 += 1
            elseif global_index_σ2 < length(σ02)
                global_index_σ1 = 1
                global_index_γ3 = 1
                global_index_γ2 = 1
                global_index_γ1 = 1
                global_index_Γ2 = 1
                global_index_Γ1 = 1
                global_index_Δ2 = 1
                global_index_Δ1 = 1
                global_index_Ω2 = 1
                global_index_Ω1 = 1
                global_index_t2 = 1
                global_index_t1 = 1
                global_index_σ2 += 1
            else
                println("All combinations exhausted. Terminating integrator.")
                finish=true
                return
            end
        
            integrator.p[1][1] = Ω01[global_index_Ω1]
            integrator.p[1][2] = Ω02[global_index_Ω2]
            integrator.p[2][1] = Δ01[global_index_Δ1]
            integrator.p[2][2] = Δ02[global_index_Δ2]
            integrator.p[3][1] = σ01[global_index_σ1]
            integrator.p[3][2] = σ02[global_index_σ2]
            integrator.p[4][1] = t00[global_index_t1]
            integrator.p[4][2] = t01[global_index_t2]
            integrator.p[5][1] = Γ01[global_index_Γ1]
            integrator.p[5][2] = Γ02[global_index_Γ2]
            integrator.p[6][1] = γ01[global_index_γ1]
            integrator.p[6][2] = γ02[global_index_γ2]
            integrator.p[6][3] = γ03[global_index_γ3]


            reinit!(integrator, integrator.u)
        end

event_cb = DiscreteCallback(condition, affect!)

The system that needs to be solved is written here

p = (Ω, Δ, σ, t0, Γ, γ)

function LambdaSystem(dρ, ρ, p, t)
    Ω₀, Δ, σ, t₀, Γ, γ = p
  
    Ω_ip = zeros(2, 1)
    
    Ω_ip = Ω₀' .* exp.(-((t .- t₀) ./ σ).^2) # Rabi Frequency
  
    Ω_ip[1] = Ω₀[1] * exp(-((t - t₀[1]) / σ[1])^2) 
    Ω_ip[2] = Ω₀[2] * exp(-((t - t₀[2]) / σ[2])^2)
    
    dρ[1] = Γ[1] * ρ[5] - Ω_ip[1] * 1im * (ρ[4] - ρ[2]) 
    dρ[5] = -(Γ[1] + Γ[2]) * ρ[5] - 1im * (Ω_ip[1] * (ρ[2] - ρ[4]) + Ω_ip[2] * (ρ[8] - ρ[6])) 
    dρ[9] = Γ[2] * ρ[5] - 1im * Ω_ip[2] * (ρ[6] - ρ[8])

    dρ[2] = (-γ[1] + 1im * Δ[1]) * ρ[2] - 1im * ((Ω_ip[1]) * (ρ[5] - ρ[1]) - Ω_ip[2] * ρ[3]) 
    dρ[3] = (-γ[2] + 1im * (Δ[1] - Δ[2])) * ρ[3] - 1im * (Ω_ip[1] * ρ[6] - Ω_ip[2] * ρ[2]) 
    dρ[6] = (-γ[3] - 1im * Δ[2]) * ρ[6] - 1im * (Ω_ip[2] * (ρ[9] - ρ[5]) + Ω_ip[1] * ρ[3])
    dρ[4] = conj(dρ[2])
    dρ[7] = conj(dρ[3])
    dρ[8] = conj(dρ[6])
end

This is a block that saves correct decisions and changes the value of the first parameter

while true
    global global_index_t1, global_index_t2, global_index_Ω1, global_index_Ω2, global_index_Δ1, global_index_Δ2
    global global_index_Γ1, global_index_Γ2, global_index_γ1, global_index_γ2, global_index_γ3
    global global_index_σ1, global_index_σ2, finish, successful_combinations

    u0 = [1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 
      0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]
    tspan = (0.0, 25.0)
    prob = ODEProblem(LambdaSystem, u0, tspan, p)
    event_cb = DiscreteCallback(condition, affect!)
    sol = solve(prob, Tsit5(), callback=event_cb)
    println(#"Event triggered at t = ", integrator.t,
    "sol",
    global_index_t1,
     global_index_t2,
     global_index_Ω1,
    global_index_Ω2,
     global_index_Δ1,
     global_index_Δ2,
     global_index_Γ1,
     global_index_Γ2,
     global_index_γ1,
    global_index_γ2,
     global_index_γ3,
     global_index_σ1,
     global_index_σ2)
    if finish==true
        break
    end
    push!(successful_combinations, 
    [Ω01[global_index_Ω1],
        Ω02[global_index_Ω2],
        Δ01[global_index_Δ1],
        Δ02[global_index_Δ2],
        σ01[global_index_σ1],
        σ02[global_index_σ2],
        t00[global_index_t1],
        t01[global_index_t2],
        Γ01[global_index_Γ1],
        Γ02[global_index_Γ2],
        γ01[global_index_γ1],
        γ02[global_index_γ2],
        γ03[global_index_γ3]])
        global_index_t1+=1
    
end

println("Successful combinations: ", successful_combinations)
@ChrisRackauckas
Copy link
Member

I don't understand what is supposed to be the bug in the package here? Can you simplify this to a clearer example? This has a bunch of other stuff going on. What is the ODE that is doing what thing wrong? That does not seem to be described in the issue.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants