0

This question contained the problem of NullSpace, but previous ones not.

This is a problem seems like my previous one, but there is some details different. I find that perhaps one method can not solve all the cases with different parameter. And thus, this question is not a duplicate of another question.

t is a generated matrix containg kz. kz are obtained by solve Det[t]==0. I find that when s and m are set as some cases (e.g. s=5,m=2), I substitue the NSolve[Numerator[Together[N@Det[Chop[t]]]] == 0, kz] into t, of course Det[t]==0, this means there must exist NullSpace[t]. However, NullSpace[t] gives a {},Why? How to solve this problem. Many Thanks! The codes are as following:

Clear["`*"]
s = 5;
m = 2;
th = Pi/4;
fi = Pi/6;
vh = 16;
mu = 11;
delta = 8;
HBAR = SetPrecision[1.05457266*10^(-34), 100];
ME = SetPrecision[9.1093897*10^(-31), 100];
ELEC = SetPrecision[1.60217733*10^(-19), 100];
Kh = Rationalize[0.211, 0];
vKh[1] := {0, 0, 0}
vKh[2] := {Kh, 0, 0};
vKh[3] := {-Kh, 0, 0};
vKh[4] := {0, Kh, 0};
vKh[5] := {0, -Kh, 0};
vKh[0] := {0, -Kh, 0};
vKh[i_] := vKh[Mod[i, 5]];
kc = Sqrt[2*ME*ELEC/HBAR^2]*10^(-11);
ku := kc*Sqrt[mu + delta];
kd := kc*Sqrt[mu - delta];
a3 = {Pi/Kh, Pi/Kh, Sqrt[2]*Pi/Kh};
k := {-ku*Sin[th]*Cos[fi], -ku*Sin[th]*Sin[fi], kz};
vkz[i_] := 
  If[Mod[i, 5] != 0, {0, 0, (i - Mod[i, 5])/5*Kh*Sqrt[2]/(m + 1)}, {0,
     0, (i - Mod[i, 5] - 5)/5*Kh*Sqrt[2]/(m + 1)}];
f[i_, i_] := 
  Total[(k + vKh[i])^2] - ku^2 - kz^2 + (kz + Total[vkz[i]])^2;
f[i_, j_] := 
  If[i == j, f[i, i], 
   kc^2*vh*Total[
     Table[Exp[I*n*Total[(vKh[j] + vkz[j] - vKh[i] + vkz[i])*a3]], {n,
        0, m}]]];
t := Array[f, {5*s, 5*s}];
slu := Select[
   kz /. NSolve[Numerator[Together[N@Det[Chop[t]]]] == 0, kz], 
   Re[#] >= 0 && Im[#] >= 0 &];
td[i_] := t /. kz -> slu[[i]];
nu[i_] := NullSpace[Chop[td[i]]];

Det[Chop[td[1]]]
nu[1]
user59546
  • 115
  • 1
  • 6
  • 2
    I think you're losing precision with machine numbers. Changing to higher precision seems to get things to work e.g. HBAR = SetPrecision[1.05457266*10^(-34), 100] and so forth. After we take the determinant, we can go back to machine precision and call NSolve on the numerator i.e. NSolve[Numerator[Together[N@Det[Chop[t]]]] == 0, kz]. – Greg Hurst Aug 06 '18 at 03:38
  • 3
    You have constants that differ by dozens of orders of magnitude; numerical results will be catastrophically bad. Please, work in natural units, where all parameters are of order one. – AccidentalFourierTransform Aug 06 '18 at 04:00
  • @AccidentalFourierTransform I mentioned that here but I don't think it's gonna happen and might not even be enough to help – b3m2a1 Aug 06 '18 at 04:02
  • @Chip Hurst Thank you very much! Now I meet another quesion. I have updated above. – user59546 Aug 06 '18 at 04:09
  • You may want to work with exact precision here. Or perhaps look into the Tolerance option for NullSpace. – Greg Hurst Aug 06 '18 at 04:24
  • Eigenvalues@td[1] show no eigenvalues close to zero. It seems Det[td[1]] != 0!? – Michael E2 Aug 06 '18 at 04:44
  • @MichaelE2 They are the following matrix which have nothing to do with the problem. – user59546 Aug 06 '18 at 04:46
  • @MichaelE2 Det[Chop[td[1]]] is a very small number (can be considered as 0) which can be verified by my codes. – user59546 Aug 06 '18 at 04:48
  • @ChipHurst I have tried nu[i_] := NullSpace[Chop[td[i]], Tolerance -> 10^(-30)], but nu[1] still gives a {}, why? – user59546 Aug 06 '18 at 04:56
  • 2
    Possible duplicate of The solutions of Det[t]==0 do not satisfy the equation – Henrik Schumacher Aug 06 '18 at 06:45
  • "Det[Chop[td[1]]] is a very small number (can be considered as 0)" -- I disagree. A nonzero number is not large or small except in relation to another. If Det[Chop[td[1]]] is zero, then at least one of the eigenvalues should be very small, say, less than 10^-8 times the largest eigenvalue if calculated at machine precision. The largest eigenvalue of td[1] is less than 1.07 in abs. val. and the smallest is greater than 0.014, which is not very small. The determinant is the product of 25 such numbers will range in magnitude from 1 down to 10^-50; it's roughly 10^-27 in your case. – Michael E2 Aug 06 '18 at 13:43
  • In the proposed duplicate, there is no substantial numerics issues. In this one, there is one, namely, underflow if Det[t] is performed at machine or other insufficient precision. Is that not distinctive enough? – Michael E2 Aug 06 '18 at 19:41
  • @AccidentalFourierTransform I don't see that difference in the magnitudes of the constants is an issue in this case. The nonzero entries of t turn out to all have a similar size. I think the problem is in the numerics of Det[t], which is connected to the size of the entries, which in turn depends on the constants. – Michael E2 Aug 06 '18 at 20:03
  • @MichaelE2 You are right. I should say Det[Chop[td[1]]] is a very small number compared with the elements (can be considered as 0). – user59546 Aug 07 '18 at 01:29

2 Answers2

4

First, I made some changes, including making everything exact:

HBAR = Rationalize[1.05457266*10^(-34), 0];
ME = Rationalize[9.1093897*10^(-31), 0];
ELEC = Rationalize[1.60217733*10^(-19), 0];

t = ExpToTrig@Array[f, {5*s, 5*s}];  (* ExpToTrig eliminates N::meprec warning *)

A warning arises when (1 + E^(-((2 I π)/3)) + E^((2 I π)/3)) is numerically evaluated at arbitrary precision. Symbolically it is zero, and ExpToTrig is one way to convert it to zero (and much faster than Simplify[t]).

Next, the following finds the nullspace:

wp = 100;
nt = N[t, wp];                             (* numericize matrix  t  *)
nd = Det[nt]; // AbsoluteTiming            (* symbolic det. with numeric coeffs. *)
With[{obj = Numerator@Together@nd},
   sol = NSolve[obj == 0 && 0 <= Re@kz && 0 <= Im@kz, kz, 
     WorkingPrecision -> MachinePrecision] (* or WorkingPrecision->Precision@obj *)
   ]; // AbsoluteTiming
sol = If[Precision@sol === MachinePrecision,
  Select[sol, Abs@First@Ratios@Eigenvalues[nt /. #][[{1, -1}]] < 10^-13 &],
  Select[sol, Precision[#] > 10 &]
  ]
NullSpace[nt /. #] & /@ sol
(*
  {0.409474, Null}              <-- timing of Det[nt]

  {0.132046, Null}              <-- timing of NSolve[]

  {{kz -> 0.105997 + 0.00517193 I}, {kz -> 0.414762 + 0.178031 I}} <-- sol

  {{{-0.0525609 + 0.0371921 I,..., -0.00461674 - 0.00426698 I}},   <-- null spaces
   {{-0.631073  + 0.195673 I,...,  -0.056275   + 0.0518164 I}}}    <--  "     "
*)

NSolve might (and does) produce spurious solutions. I used Select to eliminate them. If using MachinePrecision, the spurious ones won't have a (nearly) zero eigenvalue. If using arbitrary precision, the only spurious ones I found had an extremely low precision compared to Precision@obj, which is around 90.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Thank you very much! You are my idol! If there appear other question, I will seek guidance from you. – user59546 Aug 07 '18 at 01:24
  • @user59546 You're welcome. You might prefer the high-precision result (WorkingPrecision -> Precision@obj). The numerics of Det[t] seem a little suspect. But given the actual limited precision of the physical constants h-bar etc., you may have to put up with the rough approximation. – Michael E2 Aug 07 '18 at 03:06
  • I do not care the rough approximation, but if the machine precision can not deal with the computation, e.g. NSolve[Det[t]==0,kz] and NullSpace[t] give a result of {}. – user59546 Aug 07 '18 at 03:16
  • Now, I have changed all the allowable := to =, according to your guidance. But when I set s more than 10, e.g. s=15 and m=15, the computation can not give a result because Nomemory. I find that NSolve[obj == 0 && 0 <= Re@kz && 0 <= Im@kz, kz, WorkingPrecision -> MachinePrecision] will need too much time. How to solve this problem? Thanks a lot. – user59546 Aug 07 '18 at 03:23
  • I also find that when m is even number, the computation become much slower, e.g. s=8 and m=14 can not give sol, but s=8 and m=15 give sol quickly. – user59546 Aug 07 '18 at 03:38
  • If NullSpace[matrix] fails, you might try Eigensystem[matrix, -1]. It will give you the eigenvector and value corresponding to the eigenvalue closest to zero. If the eigenvalue is not small (compared to, say, Norm[matrix]), then maybe NSolve was not sufficiently accurate. – Michael E2 Aug 07 '18 at 13:43
  • Not NullSpace[matrix] fails, but sol = If[Precision@sol === MachinePrecision, Select[sol, Abs@First@Ratios@Eigenvalues[nt /. #][[{1, -1}]] < 10^-13 &], Select[sol, Precision[#] > 10 &]] can not be given. – user59546 Aug 08 '18 at 04:01
1

This is basically due to @MichaelE2, but I am making a few changes which might remove potential trouble spots.

One change is to set precision high, then rationalize. It might give the same results as rationalizing with second arg set to 0, I didn't check. If not, it might be a hair more faithful to the decimals. Probably makes no difference in the long run either way.

s = 5;
m = 2;
th = Pi/4;
fi = Pi/6;
vh = 16;
mu = 11;
delta = 8;
HBAR = Rationalize[SetPrecision[1.05457266*10^(-34), 100]];
ME = Rationalize[SetPrecision[9.1093897*10^(-31), 100]];
ELEC = Rationalize[SetPrecision[1.60217733*10^(-19), 100]];
Kh = Rationalize[0.211, 0];
vKh[1] = {0, 0, 0};
vKh[2] = {Kh, 0, 0};
vKh[3] = {-Kh, 0, 0};
vKh[4] = {0, Kh, 0};
vKh[5] = {0, -Kh, 0};
vKh[0] = {0, -Kh, 0};
vKh[i_] := vKh[Mod[i, 5]];

kc = Sqrt[2*ME*ELEC/HBAR^2]*10^(-11);
ku = kc*Sqrt[mu + delta];
kd = kc*Sqrt[mu - delta];
a3 = {Pi/Kh, Pi/Kh, Sqrt[2]*Pi/Kh};
k = {-ku*Sin[th]*Cos[fi], -ku*Sin[th]*Sin[fi], kz};
vkz[i_] := 
  If[Mod[i, 5] != 0, {0, 0, (i - Mod[i, 5])/5*Kh*Sqrt[2]/(m + 1)}, {0,
     0, (i - Mod[i, 5] - 5)/5*Kh*Sqrt[2]/(m + 1)}];
f[i_, i_] := 
  Total[(k + vKh[i])^2] - ku^2 - kz^2 + (kz + Total[vkz[i]])^2;
f[i_, j_] := 
  kc^2*vh*Total[
    Table[Exp[I*n*Total[(vKh[j] + vkz[j] - vKh[i] + vkz[i])*a3]], {n, 
      0, m}]];
t = Array[f, {5*s, 5*s}];

Instead of computing Det with the matrix numericized, we'll interpolated with the variable kz numericized. The advantage is this will avoid creation of denominators which in turn allows us to bypass doing polynomial algebra (e.g. Together) on an approximate rational function. I use Quiet due to the hidden zeros issue already noted. One could instead use the ExpToTrig tactic. In a production code setting I would interpolate at points on the unit circle. The below quick-and-dirty interpolation at integers is perhaps easier to follow.

Quiet[dets = Table[Det[N[t, 200]], {kz, -25, 25}]];
det = InterpolatingPolynomial[Transpose[{Range[-25, 25], dets}], kz];

We now can extract the roots and select the ones of interest. Probably better to do at high precision but I wanted to show that machine arithmetic suffices for purposes of next computing null vectors.

ns = 
 Select[NSolve[det, kz, WorkingPrecision -> MachinePrecision], 
  With[{val = kz /. #}, Re[val] >= 0 && Im[val] >= 0] &]

(* Out[390]= {{kz -> 0.105997 + 0.00517193 I}, {kz -> 
   0.414762 + 0.178031 I}} *)

td[i_] := t /. ns[[i]];

NullSpace[td[1]]

(* Out[300]= {{-0.0491714 - 0.0415701 I, -0.00914742 + 
   0.000528511 I, -0.0373394 + 0.00133131 I, -0.0108843 + 
   0.000614003 I, -0.022615 + 0.00106742 I, -0.391521 + 0.0582159 I, 
  0.0118632 - 0.420233 I, 0.00762543 - 0.0142246 I, 
  0.028503 - 0.0612858 I, 
  0.00949882 - 0.017915 I, -0.0176812 + 0.0501965 I, -0.0550383 + 
   0.204854 I, -0.020311 + 0.0653341 I, -0.0406274 + 
   0.141211 I, -0.0234697 + 0.0763275 I, 
  0.631418 + 0.247233 I, -0.0681042 - 0.324485 I, 
  0.0121311 - 0.00123308 I, 0.0563468 - 0.0137243 I, 
  0.0153808 - 0.00171599 I, -0.0300175 + 0.00143636 I, 
  0.00417262 - 0.00790418 I, 0.002762 - 0.00512705 I, 
  0.0037664 - 0.00709248 I, 0.00297463 - 0.00553832 I}} *)
Daniel Lichtblau
  • 58,970
  • 2
  • 101
  • 199