I need to solve the following equation: $$ \begin{pmatrix} \frac{\omega^2}{c^2}\varepsilon_x-\mu_z^{-1}k_y^2-\mu_y^{-1}k_z^2 & \mu_z^{-1}k_xk_y & \mu_y^{-1}k_xk_z\\ \mu_z^{-1}k_xk_y &\frac{\omega^2}{c^2}\varepsilon_y-\mu_z^{-1}k_x^2-\mu_x^{-1}k_z^2 &\mu_x^{-1}k_yk_z\\ \mu_y^{-1}k_xk_z & \mu_x^{-1}k_yk_z&\frac{\omega^2}{c^2}\varepsilon_z-\mu_y^{-1}k_x^2-\mu_x^{-1}k_y^2 \end{pmatrix} \begin{pmatrix}E_x\\E_y\\E_z\end{pmatrix}=\begin{pmatrix}0\\0\\0\end{pmatrix} $$

All parameters are given except $k_z$, $E_x$, $E_y$, $E_z$. I first find $k_z$ that makes the determinant of the matrix vanishes, then I find vector $E$ as the corresponding null space of the matrix. I have the following codes that generate the matrix and find the values of $k_z$:

import scipy as sp   # np.sqrt doesn't handle complex number?

def maxwell_matrix(k0,kx,ky,kz,eps,mu):
    kx2,ky2,kz2,k02 = kx**2,ky**2,kz**2,k0**2
    mat = sp.empty((3,3),dtype=complex)
    mat[0,0] = k02*eps[0] - ky2/mu[2] - kz2/mu[1]
    mat[1,1] = k02*eps[1] - kz2/mu[0] - kx2/mu[2]
    mat[2,2] = k02*eps[2] - kx2/mu[1] - ky2/mu[0]
    mat[0,1] = kx*ky/mu[2]
    mat[1,0] = mat[0,1]
    mat[0,2] = kx*kz/mu[1]
    mat[2,0] = mat[0,2]
    mat[1,2] = kz*ky/mu[0]
    mat[2,1] = mat[1,2]
    return mat

def compute_kz(k0,kx,ky,eps,mu):
    eps = sp.array(eps)
    mu = sp.array(mu)
    if len(eps.shape)==1:
        eps = eps[sp.newaxis,:]
    elif len(eps.shape)==0:
        eps = eps*sp.ones((1,3))
    if len(mu.shape)==1:
        mu = mu[sp.newaxis,:]
    elif len(mu.shape)==0:
        mu = mu*sp.ones((1,3))
    kx2, ky2, k02 = kx**2, ky**2, k0**2

    A = kx2*(eps[:,2]*mu[:,0]+eps[:,0]*mu[:,2])\
    B = (kx2*mu[:,0] + mu[:,1]*(ky2-k02*eps[:,2]*mu[:,0]))\
        *(kx2*eps[:,0] + eps[:,1]*(ky2-k02*eps[:,0]*mu[:,2]))

    epsz_muz = eps[:,2]*mu[:,2]
    C = sp.sqrt(A**2 - 4*epsz_muz*B)

kz1 = sp.sqrt(0.5*(-A+C)/epsz_muz)
kz2 = sp.sqrt(0.5*(-A-C)/epsz_muz)

return sp.transpose(sp.array([kz1,kz2,-kz1,-kz2]))

The equations for calculating $k_z$ are found using Mathematica and it should be correct as I have tested it many times.

I have tried two methods for finding the null space. The first method is SVD decomposition, and the second one is to find the eigenvector with eigenvalue zero. The following code does this:

import numpy as np
import scipy.linalg as spl

k0 = 1.0
theta = 37*sp.pi/180  # just some random number
phi = 0*sp.pi/180
eps = (2*np.random.rand(3)+1) + 0.1j*np.random.rand(3)
mu = (np.random.rand(3)+1)

kx = k0*sp.sin(theta)*sp.cos(phi)
ky = k0*sp.sin(theta)*sp.sin(phi)

kz = compute_kz(k0,kx,ky,eps,mu)[0]

tol = 1e-9
for i in range(4):
    mat = maxwell_matrix(k0,kx,ky,kz[i],eps,mu)  
    print('det:',spl.det(mat))   # this confirms kz is correct

    # SVD decomposition
    U, s, Vh = spl.svd(mat)
    idx = sp.argmin(sp.absolute(s))
    if sp.absolute(s[idx]) < tol:
        evect = Vh[idx]
        evect = evect/spl.norm(evect)
        raise ValueError("no solution found!")
    dot_prod = sp.dot(mat,evect)
    print('SVD vector:', evect)
    print('SVD dot prod: ({x[0]:.3e}  {x[1]:.3e}  {x[2]:.3e})    abs:{y:.3e}'.format(x=dot_prod, y=spl.norm(dot_prod)))

    # find eigenvector with zero eigenvalue
    w,v = sp.linalg.eig(mat)
    idx = sp.argmin(sp.absolute(w))
    if sp.absolute(w[idx]) < tol:
        evect = v[:,idx]
        evect = evect/spl.norm(evect)
        raise ValueError("no solution found!")
    dot_prod = sp.dot(mat,evect)
    print('method2 vector:', evect)
    print('method2 dot prod: ({x[0]:.3e}  {x[1]:.3e}  {x[2]:.3e})    abs:{y:.3e}'.format(x=dot_prod, y=spl.norm(dot_prod)))

The strange thing is that when $\varepsilon$ is real, both methods seem to give the same answer, but when $\varepsilon$ is complex, then SVD decomposition seems to fail. It agrees with method 2 for one set of $k_z$, but for the other set of $k_z$ the computed dot product is clearly zero. What's wrong with SVD decomposition/my method for finding null space? What's the best way to find the null space of a $3\times 3$ matrix, which is in general complex?

  • 2
    $\begingroup$ Please note that the SVD gives you the Hermitian transpose of V instead of V itself. Try transposing and complex conjugation first. Also please note that the singular values are ordered, thus you can leave out the argmin call. $\endgroup$ – AlexE Sep 30 '18 at 19:54

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.