Использование сложных чисел в ODE Проблема возвращает неточную ошибку

Question

Я пытаюсь реализовать уравнение Swing для системы n-Machine, используя Джулию. Когда я запускаю следующий код, я получаю это сообщение об ошибке:

LoadError: InexactError: Float64(0.0 + 1.0im)
in expression starting at /home/anna/Documents/MASTER_ARBEIT/self synchronization/first_try.jl:61
Swing_Equation(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64) at complex.jl:37
ODEFunction at diffeqfunction.jl:219 [inlined]
initialize!

Проблема возникает, поскольку я использую du[3] = (u[3] * u[2]) * im, который не может быть Float64 типом. Код работает нормально, когда я удаляю im - но тогда это уже не та модель, которую я хочу реализовать.

Как можно обойти мою проблему?

using Plots
using DifferentialEquations
inspectdr()

# Constants
P_m0  = 0.3                      # constant Mechanical Power
P_emax = 1
H = 1.01                         # Inertia constant of the system
θ_0 = asin(P_m0 / P_emax)        # angle of the system
ω_0 = 1.0                        # initial angular velocity
M = 2 * H / ω_0
D = 0.9                          # Damping constant

u02 = [θ_0;ω_0]                  # Initial Conditions
tspan = (0.0,100.0)              # Time span to solve for
p = [M;P_m0;D]  
i = 3

function Swing_Equation(du,u,t,p)                # u[1] = angle θ
    du[1] = u[2]                                 # u[2] = angular velocity ω
    P_e = real(u[3] * conj(i))
    du[2] = (1 / M) * ( P_m0 - P_e - D * u[2])   # du[2] = angular acceleration
    du[3] = (u[3] * u[2]) * im
end

# solving the differential equations

prob2 = ODEProblem(Swing_Equation,u0,tspan,p)
print(prob2)
sol2 = solve(prob2)

# Visualizing the solutoins
plot(sol2; vars = 1, label = "Θ_kura", line = ("red"))
plot!(sol2; vars = 2, label = "ω_kura", line = ("blue"))
gui()

plot(sol2,vars = (1,2),label="Kurmamoto" ,line = ("purple"))
xlabel!("Θ")
ylabel!("ω")
gui()

Answer 1

Проблема, скорее всего, в вашем вводе.

prob2 = ODEProblem(Swing_Equation,u0,tspan,p)

Я предполагаю, что в этой части вы предоставляете массив Float64 для u0? Ваш Swing_Equation затем получает u как тип Array{Float64}. Я подозреваю, что это также означает, что du то же самое.

Это приводит к сбою выражения

du[3] = (u[3] * u[2]) * im

, потому что вы пытаетесь присвоить Complex{Float64} число du[3], которое типа Float64. Затем Джулия попытается выполнить

convert(Float64, (u[3] * u[2]) * im)

, что приведет к неточной ошибке, поскольку вы не можете преобразовать комплексное число в число с плавающей запятой.

Решение состоит в том, чтобы убедиться, что du и u являются комплексными числами, поэтому вы избегаете этого преобразования. Быстрый и грязный способ решения этой проблемы - написать:

prob2 = ODEProblem(Swing_Equation, collect(Complex{Float64}, u0),tspan,p)

Это соберет все элементы в u0 и создаст новый массив, каждый элемент которого будет Complex{Float64}. Однако это предполагает 1D массив. Я не знаю твоего дела. Я не работаю с ODE-решателями самостоятельно.

Общие советы, чтобы избежать подобных проблем

Добавьте еще несколько утверждений типа в ваш код, чтобы убедиться, что вы получаете ожидаемые данные. , Это поможет выявить проблемы такого рода и упростит вам понимание происходящего.

function Swing_Equation(du::AbstractArray{T}, u::AbstractArray{T}, t,p) where T<:Complex               # u[1] = angle θ
    du[1] = u[2] :: Complex{Float64}
    P_e = real(u[3] * conj(i))
    du[2] = (1 / M) * ( P_m0 - P_e - D * u[2])   # du[2] = angular acceleration
    du[3] = (u[3] * u[2]) * im
end

Имейте в виду, что Джулия немного сложнее в сопоставлении типов, чем другие динамики c языки. Это то, что дает ему производительность.

Почему Джулия отличается от Python в этом случае?

Джулия не обновляет типы, такие как Python, до какой-либо степени. Массивы набраны. Они не могут содержать ничего, как в Python и других динамических c языках. Например, если вы создали массив, в котором каждый элемент является целым числом, вы не можете присвоить значения с плавающей запятой каждому элементу без предварительного явного преобразования в число с плавающей запятой. В противном случае Джулия должна предупредить вас, что вы получаете неточную ошибку, выдавая исключение.

В Python это не проблема, поскольку каждый элемент в массиве может иметь различный тип. Если вы хотите, чтобы каждый элемент в массиве Julia имел другой тип чисел, вы должны создать массив как Array{Number}, но это очень неэффективно.

Надеюсь, это поможет!