# plotting – 1D Integral always evaluates to 0 and Replace All

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 2×2 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.

Posted on Categories Articles