2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 Решение в Julia с компл. коэфф. Фурье: не работает nlsolve
Сообщение04.04.2025, 18:57 
Аватара пользователя
Здравствуйте!

Я решаю систему алгебраических уравнений, полученную из модели $(θ(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


 
 
 [ 1 сообщение ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group