HomotopyContinuation: Duffing Revisited: Plotting Solution from Multiple Paths

[Domenico_Lahaye]

I am seeking to plot the solutions that HomotopyContinuation.jl returns in case of multiple solution branches.

Same question reformulated: I am seeking to reproduce BifurcationKit: Duffing Equation: Time-Harmonic: Single Harmonic with HomotopyContinuation.jl

Below are two MWEs. The first is linear and works fine. The second is non-linear and fails (likely to my erroneous post-processing). The correct solution is provided by HarmonicBalance.jl by redefining the plot function (details sadly missing).

First MWE

1/ Define linear and non-linear system

@var A B p

m=1; freq = .5; om0=2*π*freq; ga=0.1; stiffnlin=1e-3; F0 = 1.;

Flin = System(
    [
        m*(om0^2-p[1]^2)*A + ga*p[1]*B,
        -ga*p[1]*A + m*(om0^2-p[1]^2)*B - F0,        
    ];
    parameters = [p]
)

F = System(
    [
       m*(om0^2-p[1]^2)*A + ga*p[1]*B + stiffnlin*.75*(A^3 + A*B^2),
       -ga*p[1]*A + m*(om0^2-p[1]^2)*B + stiffnlin*.75*(B^3 + A^2*B) - F0,
        
    ];
    parameters = [p]
)

2/ Solve and plot the linear system

# solve system for multiple values of driving frequency using start_parameters and target_parameters
omdvec = Vector(0:0.01:10.)
Avec = zeros(length(omdvec))
Bvec = zeros(length(omdvec))
for (i,p) in enumerate(omdvec) 
    # Compute all solutions for F_p
    results = solve(Flin, target_parameters = [p])
    Avec[i] = real_solutions(results)[1][1]
    Bvec[i] = real_solutions(results)[1][2]
end

plot(omdvec,Avec)
plot!(omdvec,Bvec)
plot!(omdvec,sqrt.(Avec.^2+Bvec.^2))

Second MWE

3/ Solve and plot the non-linear system

Replace “Flin” by “F” in solve function above.

PS (for what it is worth): The page How to solve a system of polynomial equations contains a tiny typo as realsolutions (without hyphen) should be real_solutions (with hyphen) instead.

[Domenico_Lahaye]

I did improved using

# solve system for multiple values of driving frequency using start_parameters and target_parameters
omdvec = Vector(0:0.01:10.)
Avec = Float64[]
Bvec = Float64[] 
start_solution = zeros(2)
start_p = 0. 
for (i,p) in enumerate(omdvec) 
    # Compute all solutions for F_p
    results = solve(F, start_solution; target_parameters = [p])
    if !(isempty(real_solutions(results))) 
        push!(Avec, real_solutions(results)[1][1]) 
        push!(Bvec, real_solutions(results)[1][2]) 
        start_solution = [real_solutions(results)[1][1], real_solutions(results)[1][2]]
        start_p = p  
    end 
end

plot(Avec)
plot!(Bvec)
plot!(sqrt.(Avec.^2+Bvec.^2))

to obtain

I am, however, not quite there yet.

[Domenico_Lahaye]

I miss functionality to plot multiple solutions for particular parameter values.

A solution is provided in HarmonicBalance.jl at HarmonicBalance.jl/src/plotting_Plots.jl at master · NonlinearOscillations/HarmonicBalance.jl · GitHub

Are there other solutions? Thx.

[Domenico_Lahaye]

Using

Avec = Vector{Float64}[]
Bvec = Vector{Float64}[]
Cvec = Vector{Float64}[]
for (i,ri) in enumerate(myrealsols)
    Aveci = []
    Bveci = [] 
    Cveci = [] 
    for (j,rj) in enumerate(ri)
        Avecj = rj[1]
        Bvecj = rj[2] 
        Cvecj = sqrt(Avecj^2+Bvecj^2)
        push!(Aveci,Avecj) 
        push!(Bveci,Bvecj) 
        push!(Cveci,Cvecj) 
    end 
    push!(Avec,Aveci)
    push!(Bvec,Bveci)
    push!(Cvec,Cveci)
end 

X = [ i*ones(length(res)) for (i,res) in pairs(myrealsols)];

scatter(X, Cvec, ms=5, legend =:false)

I now obtain

matching well with the previous solution.