search

bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
302 stars 36 forks source link

Issue with solving ODEs which contain complex variables #165

Closed RuqiShi12 closed 1 month ago

RuqiShi12 commented 1 month ago

Hi, I want to use bifurcationKit for my optical system, which contains complex variables in the equations. My code:

using BifurcationKit, Parameters, Plots

function Kerr(A, p)
    @unpack δω, P0, tau, phi, P = p
    delta_omega = δω - abs(A)^2 / (P0 * tau^2)
    d = im * exp(im * phi / 2) / sqrt(tau)
    t1 = (delta_omega * im - 1 / tau) * A
    t2 = d * P
    return t1 + t2
end

par = (δω = -200 * 1e6, P0=0.5e-3, tau=1.5e-9, phi=0., P = 0.0)

A0 = 0.0 + 0.0im  # Ensure initial condition is a floating-point complex number

prob = BifurcationProblem(Kerr, A0, par, (@lens _.P); record_from_solution = (x, p) -> (A = x[1]))

br = continuation(prob, PALC(), ContinuationPar(p_min = 0.0, p_max = sqrt(10e-3)))

plot(br)

And I got an error:

ERROR: MethodError: no method matching jacobian(::BifurcationKit.var"#7#23"{@NamedTuple{…}}, ::ComplexF64)

Closest candidates are:
  jacobian(::Any, ::Real)
   @ ForwardDiff C:\Users\rushi\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:85
  jacobian(::Any, ::AbstractArray, ::AbstractArray, ::ForwardDiff.JacobianConfig{T}, ::Val{CHK}) where {T, CHK}
   @ ForwardDiff C:\Users\rushi\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:35
  jacobian(::Any, ::AbstractArray, ::AbstractArray, ::ForwardDiff.JacobianConfig{T}) where T
   @ ForwardDiff C:\Users\rushi\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:35
  ...

Stacktrace:
  [1] (::BifurcationKit.var"#6#22")(x::ComplexF64, p::@NamedTuple{…})
    @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Problems.jl:196
  [2] jacobian(pb::BifFunction{…}, x::ComplexF64, p::@NamedTuple{…})
    @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Problems.jl:72
  [3] jacobian(pb::BifurcationProblem{…}, x::ComplexF64, p::@NamedTuple{…})
    @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Problems.jl:234
  [4] _newton(prob::BifurcationProblem{…}, x0::ComplexF64, p0::@NamedTuple{…}, options::NewtonPar{…}; normN::typeof(LinearAlgebra.norm), callback::typeof(BifurcationKit.cb_default), kwargs::@Kwargs{…})
    @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Newton.jl:88
  [5] _newton
    @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Newton.jl:62 [inlined]
  [6] newton
    @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Newton.jl:138 [inlined]
  [7] iterate(it::ContIterable{…}; _verbosity::Int64)
    @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:280
  [8] iterate
    @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:250 [inlined]
  [9] continuation
    @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:474 [inlined]
 [10] continuation(prob::BifurcationProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear_algo::Nothing, bothside::Bool, kwargs::@Kwargs{})    
    @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:549
 [11] continuation(prob::BifurcationProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…})
    @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:521
 [12] top-level scope
    @ c:\Users\rushi\OneDrive - UGent\Master\Bureaublad\Photonic Ising machine\Codings\IM on chip\dynamics\Kerr_ring.jl:20
Some type information was truncated. Use `show(err)` to see complete types.`

Could you help me resolve it?

rveltz commented 1 month ago

Hi

You have to use real values for the state. This is mainly because of ForwardDiff. A simple wrapper like this

function kerr(x::AbstractVector{<:Real}, p)
    out = Kerr(complex(x[1], x[2]), p)
    return [real(out), imag(out)]
end

should do the job.

RuqiShi12 commented 1 month ago

Hi You have to use real values for the state. This is mainly because of ForwardDiff. A simple wrapper like this function kerr(x::AbstractVector{<:Real}, p) out = Kerr(complex(x[1], x[2]), p) return [real(out), imag(out)] end should do the job.

Thanks a lot! The 'ForwardDiff' error has been resolved. however, a new error occurs when I modify code to:

# solve the steady state for a Kerr nonlinear ring

using BifurcationKit, Parameters, Plots

function Kerr(A, p)
    @unpack δω, P0, tau, phi, P = p
    delta_omega = δω - abs(A)^2 / (P0 * tau^2)
    d = im * exp(im * phi / 2) / sqrt(tau)
    t1 = (delta_omega * im - 1 / tau) * A
    t2 = d * P
    return t1 + t2
end

function kerr(x::AbstractVector{<:Real}, p)
    # Convert real parts into a complex number
    complex_state = complex(x[1], x[2])

    out = Kerr(complex_state, p)  

    return [real(out), imag(out)]
end

par = (δω = -200 * 1e6, P0=0.5e-3, tau=1.5e-9, phi=0., P = 0.0)

x0 = [0.0, 0.0]  # Ensure initial condition is a floating-point complex number

prob = BifurcationProblem(kerr, x0, par, (@lens _.P); record_from_solution = (x, p)-> complex(x[1], x[2]))

br = continuation(prob, PALC(), ContinuationPar(p_min = 0.0, p_max = sqrt(10e-3)))

plot(br)

The error shows: ERROR: ArgumentError: matrix contains Infs or NaNs Stacktrace: [1] chkfinite @ C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\lapack.jl:86 [inlined] [2] getrf!(A::Matrix{Float64}; check::Bool) @ LinearAlgebra.LAPACK C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\lapack.jl:559 [3] getrf! @ C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\lapack.jl:557 [inlined] [4] #lu!#158 @ C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\lu.jl:82 [inlined] [5] lu! @ C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\lu.jl:81 [inlined] [6] #lu#164 @ C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\lu.jl:300 [inlined] [7] lu (repeats 2 times) @ C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\lu.jl:299 [inlined] [8] \(A::Matrix{Float64}, B::Vector{Float64}) @ LinearAlgebra C:\Users\rushi\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\generic.jl:1124 [9] #_#112 @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\LinearSolver.jl:83 [inlined] [10] (::DefaultLS)(J::Matrix{Float64}, rhs::Vector{Float64}) @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\LinearSolver.jl:82 [11] _newton(prob::BifurcationProblem{…}, x0::Vector{…}, p0::@NamedTuple{…}, options::NewtonPar{…}; normN::typeof(LinearAlgebra.norm), callback::typeof(BifurcationKit.cb_default), kwargs::@Kwargs{…}) @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Newton.jl:89 [12] _newton @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Newton.jl:62 [inlined] [13] newton @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Newton.jl:138 [inlined] [14] iterate(it::ContIterable{…}; _verbosity::Int64) @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:280 [15] iterate @ C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:250 [inlined] [16] continuation(it::ContIterable{…}) @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:474 [17] continuation(prob::BifurcationProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear_algo::Nothing, bothside::Bool, kwargs::@Kwargs{}) @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:549 [18] continuation(prob::BifurcationProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}) @ BifurcationKit C:\Users\rushi\.julia\packages\BifurcationKit\5kgAI\src\Continuation.jl:521 [19] top-level scope @ c:\Users\rushi\OneDrive - UGent\Master\Bureaublad\Photonic Ising machine\Codings\IM on chip\dynamics\Kerr_ring.jl:31 Some type information was truncated. Use `show(err)` to see complete types.

The error seems to be caused by an unstable numerical computation. However, I have checked my parameter values and I can get a reasonable result from another ODE solver.

Could you help me deal with it? Thanks again!

rveltz commented 1 month ago

I find

julia> BifurcationKit.jacobian(prob, x0, par)
2×2 Matrix{Float64}:
 NaN  NaN
 NaN  NaN

It seems abs(A)^2 is not differentiable. It works with

function Kerr(A, p)
    @unpack δω, P0, tau, phi, P = p
    delta_omega = δω - abs2(A) / (P0 * tau^2)
    d = im * exp(im * phi / 2) / sqrt(tau)
    t1 = (delta_omega * im - 1 / tau) * A
    t2 = d * P
    return t1 + t2
end
RuqiShi12 commented 1 month ago

I find

julia> BifurcationKit.jacobian(prob, x0, par)
2×2 Matrix{Float64}:
 NaN  NaN
 NaN  NaN

It seems abs(A)^2 is not differentiable. It works with

function Kerr(A, p)
    @unpack δω, P0, tau, phi, P = p
    delta_omega = δω - abs2(A) / (P0 * tau^2)
    d = im * exp(im * phi / 2) / sqrt(tau)
    t1 = (delta_omega * im - 1 / tau) * A
    t2 = d * P
    return t1 + t2
end

Thanks!