# Arrays of complex values with very small imaginary part
**URL:** https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621
**Category:** Numerics
**Tags:** question, linearalgebra, complex-numbers
**Created:** [April 18, 2022, 8:14am UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621 "2022-04-18T08:14:18Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@garrek](https://discourse.julialang.org/u/garrek)
#### Post date: [April 18, 2022, 8:14am UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/1 "2022-04-18T08:14:18Z")
I am performing calculations with arrays with elements that can be complex (depending on the inputs). Later in my program, I need to check if the eigenvalues of the matrix are complex or real. (real and imaginary components have a physical meaning). The problem is that sometimes the imaginary part is very very small when it turns out the matrix is real-valued.
There are several ways to solve this problem, like initializing my arrays with zeros `zeros(ComplexF64, 4, 4)`. Or I could periodically check if there are very small numbers in my arrays, but I don’t want to accidentally throw away data. Is there a preferred / efficient way to go about this in Julia? Should I initialize my arrays differently?
This is roughly how the procedure goes:
```julia
function build_matrix(args)
M = Array{ComplexF64}(undef, 4, 4)
M[1, 1] = some calculations from the arguments
M[1, 2] = more calculations
...
return M
end
```
Then later, I have a function that uses this complex-valued matrix:
```julia
function calculate(M)
v, A = eigen(M)
if isreal(v)
perform calculations
else
calculate some other thing
end
end
```
(I’m using the LinearAlgebra package)
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [April 18, 2022, 5:53pm UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/2 "2022-04-18T17:53:31Z")
> [@garrek](#):
>
> The problem is that sometimes the imaginary part is very very small when it turns out the matrix is real-valued.
See e.g.
> [@PSA: floating-point arithmetic](https://discourse.julialang.org/t/psa-floating-point-arithmetic/8678):
>
> Sometimes people are surprised by the results of floating-point calculations such as julia\> 5/6 0.8333333333333334 # shouldn't the last digit be 3? julia\> 2.6 - 0.7 - 1.9 2.220446049250313e-16 # shouldn't the answer be 0? These are not bugs in Julia. They’re consequences of the IEEE-standard 64-bit binary representation of floating-point numbers that is burned into computer hardware, which Julia and many other languages use by default. Brief explanation You can t…
For example, if you just have some roundoff errors in your construction of the complex arrays, you could get a tiny imaginary part in the eigenvalues.
If the real-ness of eigenvalues is physically meaningful, the best approach would be to construct your matrices so that this is guaranteed. e.g. construct your matrices so that they are Hermitian (and use `Hermitian(M)`) in the cases where you expect real eigenvalues. In a physical problem, often you can do this by expressing the physics in the right way…
If this is not possible, you’ll have to resort to some kind of precision test, e.g. replace `isreal(v)` with something like `norm(imag(v)) ≤ norm(imag(v)) * eps(eltype(v)) * length(v)`. But these kinds of tests are hard to make 100% robust.
---
### Author:  [@garrek](https://discourse.julialang.org/u/garrek)
#### Post date: [April 19, 2022, 2:48am UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/3 "2022-04-19T02:48:48Z")
Thank you for the response. That gives me a few things to think about. In this particular case I’m considering propagating electromagnetic waves through lossy media, so I cannot assume the initial matrices are Hermitian. But there may be other things I can think about, like symmetry properties.
---
### Author:  [@garrek](https://discourse.julialang.org/u/garrek)
#### Post date: [April 19, 2022, 3:01am UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/4 "2022-04-19T03:01:58Z")
I just noticed that your research involves waves in photonic devices. What a coincidence! But I’m just building a simple transfer matrix implementation to interpret my experimental results. It’s not the main focus of my research and not nearly as sophisticated as what your group is doing.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [April 19, 2022, 12:11pm UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/5 "2022-04-19T12:11:58Z")
> [@garrek](#):
>
> But I’m just building a simple transfer matrix implementation to interpret my experimental results.
Note that transfer matrices can quickly become numerically unstable; it is usually better to employ [scattering matrices (S-matrices)](https://en.wikipedia.org/wiki/Redheffer_star_product#Electromagnetism) instead. Also, the relevant eigenvalues of the S matrices are typically the [nonlinear eigenfrequencies](https://en.wikipedia.org/wiki/Nonlinear_eigenproblem) where \det S^{-1}(\omega) = 0 (the resonances), which all lie in the lower-half complex plane for lossy media and open systems. Which eigenvalues are you thinking of that are purely real?
---
### Author:  [@garrek](https://discourse.julialang.org/u/garrek)
#### Post date: [April 20, 2022, 6:04am UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/6 "2022-04-20T06:04:33Z")
Now I’m aware of the instabilities, but when I started my advisor gave me Pochi Yeh’s book on transfer matrices. Not being an expert I looked into more recent literature and found some papers that claim to resolve singularities and discontinuities by sorting the modes by real / complex and by sign (forward / backward propagating). The modes are the eigenmodes of a 4x4 matrix, \Delta, with elements relating to the elements of the material dielectric tensor\[1\]. Then you solve the equation
q\_{ij}\Delta(i) = \Psi\_{ij}\Delta(i)
where q\_{ij} are the modes that are sorted into real wave vectors that propagate and complex wave vectors that are damped. The main paper that I’m following\[2\] claims to combine several approaches for a generalized transfer matrix to avoid numerical instabilities.
My Julia implementation is able to reproduce the tutorial plots from these authors, and I can produce simple Fabry-Pérot fringes (that I need for my experiment), so I think my code is correct. **However, I want to rigorously test my code and I have been struggling to write good unit tests for this sorting part so that I am sure I can identify purely real eigenvalues.** Sometimes when I was refactoring my code, I got very small imaginary values and I want to know in my tests whether these are mistakes or physical results. (Of course, strange plots tell me, but I prefer if I can write a test or an error message or something).
If scattering matrices will resolve these issues and are easier to implement, I will certainly do that. I have run across them, but haven’t tried to implement S-matrices and I’m really not familiar with how they work compared to transfer matrices.
(My goal, by the way, is to simulate absorbing material in between two planar metallic layers).
* * *
1. D. W. Berreman, J. Opt. Soc. Am. **62** , 502 (1972).
2. N. C. Passler and A. Paarmann, J. Opt. Soc. Am. B **34** , 2128 (2017).
---
### Author:  [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)
#### Post date: [April 26, 2022, 3:00am UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/7 "2022-04-26T03:00:49Z")
In case it is useful to you, Chapter 3 of the [Theory Documentation](https://github.com/simonp0420/PSSFSS.jl/blob/main/docs/TheoryDocs/theorydoc.pdf) for my PSSFSS package derives the Redheffer formula for cascading S matrices from first principles.
---
### Author:  [@garrek](https://discourse.julialang.org/u/garrek)
#### Post date: [April 27, 2022, 4:46am UTC](https://discourse.julialang.org/t/arrays-of-complex-values-with-very-small-imaginary-part/79621/8 "2022-04-27T04:46:50Z")
Thank you! I will check it out.