1
$\begingroup$

I am a bit puzzled by Mathematica's behaviour. I am trying to evaluate some 1D integrals which depend on parameters phi1,phi2 numerically (functions f1 and f2), however NIntegrate always outputs 0 whereas plotting the integrand with Plot shows that the integral is non-zero.

Here is the code in question:

t = 1;
lambda = 0.3;
hbarOmega = 4;
q = 0.5;
alpha = 4 lambda^2 t (Sin[q/2])/q;

A[k_] := 2 t (Sin[q/2]/q) Sin[k] (1/(hbarOmega + 4 t Sin[q/2] Sin[k]) + 
 1/(hbarOmega - 4 t Sin[q/2] Sin[k]));

B[k_] := 2 t (Sin[q/2]/q) Sin[k] (1/(hbarOmega + 4 t Sin[q/2] Sin[k]) - 
 1/(hbarOmega - 4 t Sin[q/2] Sin[k]));

e1[k_] := 4 t lambda^2  (Sin[q/2]/q) B[k + q/2] Sin[k + q/2];
e2[k_] := -4 t lambda^2  (Sin[q/2]/q) B[k - q/2] Sin[k - q/2];
en[k_] := -2 t Cos[k] + e1[k - q] + e2[k + q];
F[k_, phi1_, phi2_] := - alpha (A[k] phi1 + Sin[k] phi2);
mu[phi1_, phi2_] := - 2 alpha Re[phi1 Conjugate[phi2]];

hMf[k_, phi1_, phi2_] := {{en[k + q/2],F[k, phi1, phi2]}, {Conjugate[F[k, phi1, phi2]], en[k - 
q/2]}};


u[k_, phi1_, phi2_] := Transpose[Eigenvectors[hMf[k, phi1, phi2]]];
enMf[k_, phi1_, phi2_] := ConjugateTranspose[u[k, phi1, phi2]].hMf[k, phi1, phi2].u[k, phi1, 
phi2];

o[k_, phi1_, phi2_] := ConjugateTranspose[u[k, phi1, phi2]].{{0, 0}, {1, 0}}.u[k, phi1,phi2];


f1[phi1_?NumericQ,phi2_?NumericQ] := (1/(2 Pi)) NIntegrate[
Sin[k]*(o[k, phi1, phi2][[1, 1]] Boole[
     enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
   o[k, phi1, phi2][[2, 2]] Boole[
     enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]), {k, -Pi, Pi}];

f2[phi1_?NumericQ,phi2_?NumericQ] := (1/(2 Pi)) NIntegrate[A[k] (o[k, phi1, phi2][[1, 1]] 
Boole[enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
   o[k, phi1, phi2][[2, 2]] Boole[
     enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]), {k, -Pi,Pi}];

f1[1, 1]

Plot[Sin[k]*(o[k, phi1, phi2][[1, 1]] Boole[
   enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
 o[k, phi1, phi2][[2, 2]] Boole[
   enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]) /. {phi1 -> 1, phi2 -> 1}, {k, -Pi, Pi}]

Plot[Sin[k]*(o[k, 1, 1][[1, 1]] Boole[enMf[k, 1, 1][[1, 1]] <= mu[1, 1]] + o[k, 1, 1][[2, 2]] 
Boole[enMf[k, 1, 1][[2, 2]] <= mu[1, 1]]), {k, -Pi, Pi}]

So for example, when I evaluate the integral given by f1[1,1], I get an output of 0 alongside with a warning that Mathematica couldn't find all eigenvectors of a 2x2 matrix which is bizarre, because when I plot the eigenvalues separately with:

Plot[{enMf[k, 1, 1][[1, 1]], enMf[k, 1, 1][[2, 2]]}, {k, -Pi, Pi}]

the warning is not present.

However when I plot the integrand inside f1 with the same parameters phi1=1,phi2=1 with:

Plot[Sin[k]*(o[k, 1, 1][[1, 1]] Boole[
  enMf[k, 1, 1][[1, 1]] <= mu[1, 1]] + 
o[k, 1, 1][[2, 2]] Boole[
  enMf[k, 1, 1][[2, 2]] <= mu[1, 1]]), {k, -Pi, Pi}]

it shows that the integrand is always positive and that the integral should be non-zero.

Moreover, there is another oddity: replacing the parameters phi1,phi2 within the integrand with /. for the same values phi1=1,phi2=1 like so:

Plot[Sin[k]*(o[k, phi1, phi2][[1, 1]] Boole[
   enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
 o[k, phi1, phi2][[2, 2]] Boole[
   enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]) /. {phi1 -> 1,phi2 -> 1}, {k, -Pi, Pi}]

gives a different plot compared to just substituting the values phi1=1,phi2=1 directly:

Plot[Sin[k]*(o[k, 1, 1][[1, 1]] Boole[
  enMf[k, 1, 1][[1, 1]] <= mu[1, 1]] + 
o[k, 1, 1][[2, 2]] Boole[
  enMf[k, 1, 1][[2, 2]] <= mu[1, 1]]), {k, -Pi, Pi}]

which I find very bizarre and assume has to do with delayed evaluation. Any suggestions would be appreciated.

$\endgroup$
3
$\begingroup$

You're right that it's a delayed evaluation problem! I traced the problem back to Eigenvectors. Note that it behaves differently when given symbolic arguments and numeric ones:

Eigenvectors[hMf[0.01, phi1, phi2]] /. {phi1 -> 1, phi2 -> 1}

(* Out: {{-4.064089631, 1.}, {0.2460575654, 1.}} *)

Eigenvectors[hMf[0.01, 1, 1]]

(* Out: {{0.2389308952, 0.9710365736}, {-0.9710365736, 0.2389308952}} *)

The deal is that Mathematica both 1) doesn't normalize eigenvectors when symbolic and 2) might order them differently. Note that

Normalize /@ (Eigenvectors[hMf[0.01, phi1, phi2]] /. {phi1 -> 1, phi2 -> 1})

(* Out: {{-0.9710365736, 0.2389308952}, {0.2389308952, 0.9710365736}} *)

Which is the same as Eigenvectors[hMf[0.01, 1, 1]], but in reverse order.

Sort[Normalize /@ (Eigenvectors[hMf[0.01, phi1, phi2]] /. {phi1 -> 1, phi2 -> 1})]
  == Sort[Eigenvectors[hMf[0.01, 1, 1]]]

(* Out: True *)

So, you'll need to decide whether the vectors should be normalized, and whether the order matters!

But there's also another issue: different normalizations can be chosen. Normalized eigenvectors are unique only up to norm 1 elements of your field, which in this case are $\pm 1$ if we're doing this over the reals. Which normalization do you want? For this reason, even after redefining u[k_, phi1_, phi2_] := Transpose[Normalize /@ Eigenvectors[hMf[k, phi1, phi2]]];, we still have

o[0.01, phi1, phi2][[2, 2]] /. {phi1 -> 1, phi2 -> 1}

(* Out: 0.2320106378 *)

o[0.01, 1, 1][[2, 2]] 
(* Out: -0.2320106378 *)

So that's another choice to make as well.


However, there is a way to avoid some of this, since I noticed you seem to be using the matrix of eigenvectors for diagonalization. Check out JordanDecomposition and other *Decomposition functions, which return not only the decomposition, but similarity matrices you can use as well! Also look at Eigensystem, which returns a list {evalues, evectors} for which M.Transpose[evectors] == Transpose[evectors].DiagonalMatrix[evalues]. You could use this to save on some computation in calculating enMf by simply using DiagonalMatrix[evalues].

$\endgroup$
3
  • $\begingroup$ Thanks for this. Concerning your remark on normalization (last bit of code in your answer), I think that its the ordering of the vectors that's causing the sign change, as the structure of the matrix $o$ ( $o = u^{\dagger} \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} u $) implies that its diagonal elements are invariant under the change in sign of any of the eigenvectors $v_1,v_2$ making up $u = (v_1,v_2)$. Also, now I understand why the plots of the integrand are different but I am still not sure why the integral f[1,1] is zero since a plot of the integrand suggests otherwise. $\endgroup$ Sep 7 '21 at 8:48
  • $\begingroup$ Update: so making all the functions from hMf down numerical with ?NumericQ gives a non-zero integral, however now I get error messages like so: "Part specification o[k,2,2][[1,1]] is longer than depth of object" which I guess arise from Part being evaluated prematurely. I am a bit confused by this message because since we are using NIntegrate, shouldnt k always have a numerical value? $\endgroup$ Sep 7 '21 at 9:26
  • 1
    $\begingroup$ @RuslanMushkaev Try Indexed[o[k.2.2], {1, 1}] instead of Part. NIntegrate inspects the integrand symbolically to determine which method to choose. Good thing you didn’t have o[k, 2, 2][[1]], because that evaluates to k, which would give you strange results and no error. $\endgroup$
    – Michael E2
    Sep 7 '21 at 14:08

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.