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:
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:
function calculate(M)
v, A = eigen(M)
if isreal(v)
perform calculations
else
calculate some other thing
end
end
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.
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.
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.
Note that transfer matrices can quickly become numerically unstable; it is usually better to employ scattering matrices (S-matrices) instead. Also, the relevant eigenvalues of the S matrices are typically the nonlinear eigenfrequencies 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?
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).
D. W. Berreman, J. Opt. Soc. Am. 62 , 502 (1972). ↩︎
N. C. Passler and A. Paarmann, J. Opt. Soc. Am. B 34 , 2128 (2017). ↩︎
In case it is useful to you, Chapter 3 of the Theory Documentation for my PSSFSS package derives the Redheffer formula for cascading S matrices from first principles.