Why the null(tuu1) gives a Empty sym: 1-by-0?

tuu is a generated matrix containing kz. Through solve the equation det(tuu)==0, one can obtain the solutions kz. I select the ones satisfy real(kz0)>=0&imag(kz0)>=0. Now, I substitute the five solutions into tuu, and of course, det(tuu1) should be zero (actually it is 8*e-15 because of the float number). According to the law of linear algebra, there must exist non-trival null(tuu1). However, the program gives the result of "Empty sym: 1-by-0". Why? How to solve this problem? Many thanks! The codes are as follows.

clear all;
syms kz x y
x=0.5;
y=0;
m = 2;
vh = 4;
mu = 11;
delta = 8;
HBAR = 1.05457266e-34;
ME = 9.1093897e-31;
ELEC = 1.60217733e-19;
Kh = 2.106;
vKh = [0,0,0;Kh,0,0;-Kh,0,0;0,Kh,0;0,-Kh,0];
kc = sqrt(2*ME*ELEC/HBAR^2)*1e-10;
ku = kc*sqrt(mu+delta);
kd = kc*sqrt(mu-delta);
a3 = [pi/Kh,pi/Kh,sqrt(2)*pi/Kh];
kuu =@(x,y) [-ku*sin(x)*cos(y), -ku*sin(x)*sin(y), kz];
n=0:m;
for p=1:5;
    for q=1:5;
    tuu(p,q)=(sum((kuu(x,y) + vKh(p,:)).^2)-ku^2)*(p==q)+ kc^2*vh*sum(exp(i*n*sum((vKh(q,:)-vKh(p,:)).*a3)))*(p~=q);
    end
end
dtuu=det(tuu);
kz0=vpasolve(dtuu,kz);
kz0=double(kz0);
kzz=kz0(real(kz0)>=0&imag(kz0)>=0);
tuu1=subs(tuu,kz,kzz(1));
tuu2=subs(tuu,kz,kzz(2));
tuu3=subs(tuu,kz,kzz(3));
tuu4=subs(tuu,kz,kzz(4));
tuu5=subs(tuu,kz,kzz(5));
null(tuu1)

Answer 1

[Walter Roberson]

The round off error associated with those floating point constants leads to your tuu1 having rank 5 and so has no null space.

If you carefully work with symbolic values, then you can get rank of tuu1 to be 4. Unfortunately in that form, there seems to be something odd or wrong about the internal mupad linalg::rank routine when you feed it the indefinitely precise version of tuu1 .

The work-around is to calculate in symbolic form as long as possible, but then to null() double(tuu1)

syms kz x y
Q = @(v) sym(v, 'r');
x=Q(0.5);
y=Q(0);
m = Q(2);
vh = Q(4);
mu = Q(11);
delta = Q(8);
HBAR = Q(1.05457266e-34);
ME = Q(9.1093897e-31);
ELEC = Q(1.60217733e-19);
Kh = Q(2.106);
vKh = [0,0,0;Kh,0,0;-Kh,0,0;0,Kh,0;0,-Kh,0];
kc = sqrt(2*ME*ELEC/HBAR^2)*Q(1e-10);
ku = kc*sqrt(mu+delta);
kd = kc*sqrt(mu-delta);
a3 = [pi/Kh,pi/Kh,sqrt(2)*pi/Kh];
kuu =@(x,y) [-ku*sin(x)*cos(y), -ku*sin(x)*sin(y), kz];
n=0:m;
tuu = zeros(5,5,'sym');
for p=1:5
    for q=1:5
    tuu(p,q)=(sum((kuu(x,y) + vKh(p,:)).^2)-ku^2)*(p==q)+ kc^2*vh*sum(exp(1i*n*sum((vKh(q,:)-vKh(p,:)).*a3)))*(p~=q);
    end
end
dtuu=det(tuu);
kz0=solve(dtuu,kz);
kz0d = double(kz0);
kzz=kz0(real(kz0d)>=0&imag(kz0d)>=0);
tuu1=subs(tuu,kz,kzz(1));
tuu2=subs(tuu,kz,kzz(2));
tuu3=subs(tuu,kz,kzz(3));
tuu4=subs(tuu,kz,kzz(4));
tuu5=subs(tuu,kz,kzz(5));
null( double(tuu1) )