Krylov-Bogoliubov for commensurate frequencies

[oameye]

Krylov-Bogoliubov method for commensurate frequencies, e.g., $\omega$ and $3\omega$. Technically the one can compute the slow-flow equations for such a system, but at the moment KB in HB.jl compute the wrong equations.

eq1 = d(d(x, t), t) + ω₁^2 * x + α₁ * x^3 + 3 * J₁ * x^2 * y + J₂ * y^2 * x #+ γ₁ * d(x, t) + η₁ * x^2 * d(x, t)
eq2 = d(d(y, t), t) + ω₂^2 * y + α₂ * y^3 + J₁ * x^3 + J₂ * x^2 * y #+ γ₂ * d(y, t) + η₂ * y^2 * d(y, t)


forces = [F * cos(ω * t), 0]

dEOM_temp = DifferentialEquation([eq1, eq2] - forces, [x, y]);

add_harmonic!(dEOM_temp, x, ω) # x will rotate at ω
add_harmonic!(dEOM_temp, y, 3 * ω); # y will rotate at 3*ω

# harmonic_tmp = get_krylov_equations(dEOM_temp, order=1);
harmonic_tmp = get_harmonic_equations(dEOM_temp);

rearranged = HB.rearrange_standard(harmonic_tmp)
prob_temp = HB.Problem(rearranged)
rearranged.equations

[oameye]

Yes KB in HB.jl does something wrong

[oameye]

        @variables t x(t) y(t) ω0 ω F α J # symbolic variables
        eq1 = d(d(x, t), t) + ω0^2 * x + α * x^3 ~ F * cos(ω * t) + J * y
        eq2 = d(d(y, t), t) + ω0^2 * y + α * y^3 ~ F * cos(ω * t) + J * x
        EOM = DifferentialEquation([eq1, eq2], [x, y])
        add_harmonic!(EOM, x, ω)
        add_harmonic!(EOM, y, 3*ω)
        krylov_eq = get_krylov_equations(EOM; order=1)
        harmonic_eq = get_harmonic_equations(EOM)
        rearranged = HarmonicBalance.rearrange_standard(harmonic_eq)


            eqk = expand_fraction(krylov_eq.equations[1].lhs)
            eqh = expand_fraction(rearranged.equations[1].lhs)
            @variables T u1(T) v1(T) ω0 ω F α # symbolic variables
            subs = Dict([u1, v1, α, F, ω0, ω, J] .=> rand(7))
            solk = substitute(eqk, subs)
            solh = substitute(eqh, subs)
            @test Float64(solk + solh) == 0.0