Complex functions

[dpsanders]

It would be great to be able to extract real and imaginary parts of a complex function.
Currently this does not seem to work, e.g.

julia> @variables a b
(a, b)

julia> c = a + b * im
a + b * (im)

julia> c^2
(a + b * (im)) ^ 2

It needs to know that im^2 is -1.

[dpsanders]

I guess this "just" needs a simplifier.

[ChrisRackauckas]

Yeah it sounds like it's a job for a simplification tool

[dpsanders]

This partially works using a custom complex type:

using ModelingToolkit

struct Compl{T}
    re::T
    im::T
end

Base.show(io::IO, x::Compl) = print(io, "($(x.re)) + i*($(x.im))")

import Base: +, -, *, real, imag, conj

+(a::Compl, b::Compl) = Compl(a.re + b.re, a.im + b.im)
*(a::Compl, b::Compl) = Compl(a.re * b.re - a.im * b.im, a.re*b.im + a.im*b.re)

conj(a::Compl) = Compl(a.re, -a.im)

@variables a b
x = Compl(a, b)

real(x::Compl) = x.re
imag(x::Compl) = x.im

The simplification is the problem:

julia> x = Compl(a, b)
(a) + i*(b)

julia> x * conj(x)
(a * a - b * -b) + i*(a * -b + b * a)


julia> y = x * conj(x)
(a * a - b * -b) + i*(a * -b + b * a)

julia> imag(y)
a * -b + b * a

cc @shashi

[dpsanders]

By the way, this is another example where Julia's restriction to defining Complex{T} only for T <: Real is incorrect.

[dpsanders]

Is this now doable with the new simplification features? @shashi

[shashi]

Your custom type works if you pass the final result to ModelingToolkit.simplify_constants (I just pushed a fix to SymbolicUtils that is needed though).

[dpsanders]

You can now simplify the real and imaginary parts separately:

julia> ex = x * conj(x)
((a * a) - (b * (-b))) + i*((a * (-b)) + (b * a))

julia> simplify(real(ex))
(a ^ 2) + (b ^ 2)

julia> simplify(imag(ex))
0

[AshtonSBradley]

so I guess it is a straightforward method extension for simplify to just add up the results of real and imag simplify?

[YingboMa]

This now works.

julia> @variables a b
(a, b)

julia> c = a + b * im
a + (b)*im

julia> c*c*c
a*((a^2) - (b^2)) - (2a*(b^2)) + (b*((a^2) - (b^2)) + 2b*(a^2))*im