Решение в Julia с компл. коэфф. Фурье: не работает nlsolve

[Rozarya]

Здравствуйте!

Я решаю систему алгебраических уравнений, полученную из модели $(θ(t))'=ω+σ \sin(θ(t))$ на основе гармонического баланса. Неизвестными являются комплексные коэффициенты Фурье Ψ_p и период T. Я написала код в Julia с использованием пакета NLsolve.

Но почему-то никакой из методов (:newton, :trust_region, :anderson) не даёт сходимости. Пробовала разные начальные приближения — результат тот же. Сначала я пыталась решить систему для системы уравнений без разделения на действительных и мнимых частей уравнений, потом разделила систему на действительные и мнимые части. Но в обеих кодах результат один тот же. Помогите пожалуйста, у меня мозг уже плавится...

Написанный мной код:

using LinearAlgebra, NLsolve, Plots, Roots, DelimitedFiles

# Differentiation
D_matrix(M, T) = Diagonal([2π * im * p / T for p in -M:M])

# Vandermonde matrix
vandermonde(M) = [exp(2π * im * p * n / (2M + 1)) for n in -M:M, p in -M:M]

# Delay matrix (Gamma)
Gamma_matrix(M, τ, T) = Diagonal([exp(-2π * im * p * τ / T) for p in -M:M])

# Coupling function
H(θ, a, r) = -sin(θ - a) + r * sin(2θ)
H_1(θ) = sin(θ)

# Construction of the Fourier matrix from x
build_Psi(x, N, M) = reshape(complex.(x[1:end-1]), 2M+1, N)

x0 = [0.1563485897289011 + 0.12169282865097303im,
 -0.27550661452108116 - 0.5640295616468156im, 
 -7.400578237139926 + 0.0im, 
 -0.27550661452108116 + 0.5640295616468156im, 
 0.1563485897289011 - 0.12169282865097303im,
14.414615569058912] 

M = 2
N = (length(x0) - 1) ÷ (2M + 1)

ω = 1.0 
σ = 0.9





# solve_system
S = vandermonde(M)
S_inv = conj(S)./ (2M + 1)

function f!(F, x)
    T = x[end]
    D = D_matrix(M, T) 
    Ψ = build_Psi(x, N, M)
    X = real(S * Ψ)
    X_shift = X .+ (2π .* collect(-M:M) ./ (2M + 1))

    lhs = S * D * S_inv * X .+ 2π / T  
    rhs = σ * H_1.((X_shift)) .+ ω
        
    
    for p in 1:(2M+1) 
        F[p] = real(lhs - rhs)[p]
    end
    
    # Average normalization condition
    F[end] = real(sum(p*Ψ[p+M+1]-2π for p in -M:M))#(sum(Ψ) - π)
    #F[end] = real(sum(Ψ[p+M+1]-2π for p in -M:M))#(sum(Ψ) - π)

    


end


solver_methods = [:newton]
tolerances = [1e-10, 1e-8, 1e-6, 1e-4, 1e-3, 1e-2, 1e-1]
#x_initial = real(x0)
for method in solver_methods
    for ftol in tolerances
        try
            res = nlsolve(f!, x0;
                method=method,
                ftol=ftol,
                xtol=ftol,
                iterations=10000)
            if res.f_converged
                println("Converged with method: $method, tolerance: $ftol")
                T = res.zero[end]
                Ψ = build_Psi(res.zero, N, M)
                
                println("T = ", T)
                println("Ψ = ", Ψ)    
                return T, Ψ 

            else
                println("Not converged")
            end
        catch e
            println("Failed with method: $method, tolerance: $ftol")
            println("Error: ", e)
        end
    end
end