Search code examples
juliascientific-computingdifferentialequations.jl

Using Complex Numbers in ODE Problem returns Inexact Error


I am trying to implement to Swing equation for a n-Machine system using Julia. When i run the following code I get this Error Message:

LoadError: InexactError: Float64(0.0 + 1.0im)
in expression starting at /home/Documents/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!

The problem is occuring since I am using du[3] = (u[3] * u[2]) * im which can not be a Float64 type. The code is working fine when I remove the im - but then it is not the model I want to implement anymore.

What way is there to work around my problem?

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()

Solution

  • The problem is most likely in your input.

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

    I am guessing that in this part you are providing an array of Float64 for u0? Your Swing_Equation then receives u as an Array{Float64} type. I suspect that also means du is the same.

    This causes the expression

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

    to fail because you are trying to assign a Complex{Float64} number to du[3] which is of type Float64. Julia will then try to perform a

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

    Which will cause the inexact error, because you cannot convert a complex number to a floating point number.

    The Solution is to make sure du and u are complex numbers so you avoid this conversion. A quick and dirty way to solve that would be to write:

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

    This will collect all elements in u0 and create a new array where every element is a Complex{Float64}. However this assumes a 1D array. I don't know your case. I don't work with ODE solvers myself.

    General Advice to avoid this kind of problem

    Add some more type assertions to in your code to make sure you get the kind of inputs you expect. This will help catch these kinds of problem and make you more easily see what is going on.

    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
    

    Keep in mind Julia is a bit more demanding when it comes to matching up types than other dynamic languages. That is what gives it the performance.

    Why is Julia different from Python in this case?

    Julia does not upgrade types like Python to whatever fits. Arrays are typed. They cannot contain anything like in Python and other dynamic languages. If you e.g. made an array where each element is an integer, then you cannot assign float values to each element without doing an explicit conversion to floating point first. Otherwise Julia has to warn you that you are getting an inexact error by throwing an exception.

    In Python this is not a problem because every element in an array can be a different type. If you want every element in a Julia array to be a different number type then you must create the array as a Array{Number} type but these are very inefficient.

    Hope that helps!