-1
$\begingroup$

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]
$\endgroup$
  • 2
    $\begingroup$ 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]. $\endgroup$ – Chip Hurst Aug 6 '18 at 3:38
  • 3
    $\begingroup$ 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. $\endgroup$ – AccidentalFourierTransform Aug 6 '18 at 4:00
  • $\begingroup$ @AccidentalFourierTransform I mentioned that here but I don't think it's gonna happen and might not even be enough to help $\endgroup$ – b3m2a1 Aug 6 '18 at 4:02
  • $\begingroup$ @Chip Hurst Thank you very much! Now I meet another quesion. I have updated above. $\endgroup$ – user59546 Aug 6 '18 at 4:09
  • 2
    $\begingroup$ Possible duplicate of The solutions of Det[t]==0 do not satisfy the equation $\endgroup$ – Henrik Schumacher Aug 6 '18 at 6:45
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you very much! You are my idol! If there appear other question, I will seek guidance from you. $\endgroup$ – user59546 Aug 7 '18 at 1:24
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Aug 7 '18 at 3:06
  • $\begingroup$ 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 {}. $\endgroup$ – user59546 Aug 7 '18 at 3:16
  • $\begingroup$ 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. $\endgroup$ – user59546 Aug 7 '18 at 3:23
  • $\begingroup$ 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. $\endgroup$ – user59546 Aug 7 '18 at 3:38
1
$\begingroup$

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}} *)
$\endgroup$
  • $\begingroup$ Thank you very much! It helps me a lot. $\endgroup$ – user59546 Aug 13 '18 at 1:43

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.