How can I perform a successive NIntegrate process with two different variables?

Question

f[x_, y_, z_] ={
  {1.5 - 0.75 (Cos[x] + Cos[y]), -I Sin[x] - Sin[y], 
   3.5 - 1.5 E^(-I z) - 1.5 (Cos[x] + Cos[y]), 0},
  {I Sin[x] - Sin[y], 1.5 - 0.75 (Cos[x] + Cos[y]), 0, 
   3.5 - 1.5 E^(-I z) - 1.5 (Cos[x] + Cos[y])},
  {3.5 - 1.5 E^(I z) - 1.5 (Cos[x] + Cos[y]), 0, 
   1.5 - 0.75 (Cos[x] + Cos[y]), I Sin[x] + Sin[y]},
  {0, 3.5 - 1.5 E^(I z) - 1.5 (Cos[x] + Cos[y]), -I Sin[x] + Sin[y], 
   1.5 - 0.75 (Cos[x] + Cos[y])}
 };    
O1[x_, y_, z_] = D[f[x, y, z], x];
O2[x_, y_, z_] = D[f[x, y, z], z];
O3 = KroneckerProduct[PauliMatrix[0], PauliMatrix[2]];
O4[x_, y_, z_] = 1/4 (O2[x, y, z].O3 + O3.O2[x, y, z]);    
Nr[fi_, r_, br_] := (Inverse[(r + I*br)*IdentityMatrix[4] - fi]);
Na[fi_, r_, br_] := (Inverse[(r - I*br)*IdentityMatrix[4] - fi]);
Fxyz[x_, y_, z_, r_, br_] := 
  Re[Tr[O1[x, y, z].Nr[f[x, y, z], r, 
       br].(-O4[x, y, z].Nr[f[x, y, z], r, br] + 
        Nr[f[x, y, z], r, br].O4[x, y, z]).Nr[f[x, y, z], r, br]]];

I have the function Fxyz(x,y,z,r,br)as defined above. I would like to perform the following process:
$$text{G}(r)=int _{-pi }^{pi }dx dy dz  text{Fxyz}(x,y,z,r,0.01)$$ and then
$$text{F}(v)=int _{-2}^vdrtext{G}(r)$$ with $v$ runs from $-1$ to $2$, at the end I would like to get a plot for $text{F}(v)$ as a function of $v$. How Can I do that?
Here is my try using NIntegrate
G[r_] := NIntegrate[
  Fxyz[x, y, z, r, 
   0.01], {x, -[Pi], [Pi]}, {y, -[Pi], [Pi]}, {z, -[Pi], [Pi]}, 
  Method -> {Automatic, "SymbolicProcessing" -> 0}]
Fr = ParallelTable[
  NIntegrate[G[r], {r, -2, v}, 
   Method -> {Automatic, "SymbolicProcessing" -> 0}], {v, -1, 2, 
   0.05}]

but the integrals are evaluated to non-numeric values?!

Hausdorff · Answer

The integrand can be simplified to a significantly more compact form.
First, the floats in f[x,y,z] can be replaced by rationals, which are easier for Mathematica to deal with:
f[x_,y_,z_] = {{3/2 - (3*(Cos[x] + Cos[y]))/4, (-I)*Sin[x] - Sin[y], 
    7/2 - 3/(2*E^(I*z)) - (3*(Cos[x] + Cos[y]))/2, 
    0}, {I*Sin[x] - Sin[y], 3/2 - (3*(Cos[x] + Cos[y]))/4, 0, 
    7/2 - 3/(2*E^(I*z)) - (3*(Cos[x] + Cos[y]))/2}, {7/
      2 - (3*E^(I*z))/2 - (3*(Cos[x] + Cos[y]))/2, 0, 
    3/2 - (3*(Cos[x] + Cos[y]))/4, I*Sin[x] + Sin[y]}, {0, 
    7/2 - (3*E^(I*z))/2 - (3*(Cos[x] + Cos[y]))/2, (-I)*Sin[x] + 
     Sin[y], 3/2 - (3*(Cos[x] + Cos[y]))/4}};

After that one can repeatedly use Simplify and FullSimplify
Nrf[x_, y_, z_, r_, br_] = Nr[f[x, y, z], r, br] // FullSimplify;

intMat[x_, y_, z_, r_, br_] = O1[x, y, z].Nrf[x, y, z, r, br].
    (Nrf[x, y, z, r, br].O4[x, y, z] - O4[x, y, z].Nrf[x, y, z, r, br]).
    Nrf[x, y, z, r, br] // Simplify // FullSimplify;

FxyzNew[x_, y_, z_, r_, br_] = 
 Tr[intMat[x, y, z, r, br]] // Simplify // Re // ExpToTrig // Simplify

which results in
FxyzNew[x,y,z,r,br]

-1536*Im[(6*Cos[z] + 2*Cos[x]*(3 + (-7 + 3*Cos[y])*Cos[z]))/
(478 + (96*I)*br + 32*br^2 + 96*r - (64*I)*br*r - 32*r^2 + 11*Cos[2*x] - 
264*Cos[y] - (48*I)*br*Cos[y] - 48*r*Cos[y] + 11*Cos[2*y] - 336*Cos[z] + 
 144*Cos[y]*Cos[z] + 12*Cos[x]*(-22 - (4*I)*br - 4*r + 9*Cos[y] + 12*Cos[z]))^2]

Comparing the timings for numerical evaluation, there is a significant improvement:
(ClearSystemCache[]; FxyzNew[1, 2, 3, 4, 5]) // N // RepeatedTiming
(ClearSystemCache[]; Fxyz[1, 2, 3, 4, 5]) // N // RepeatedTiming

{0.00060, -0.00324018}
{0.681, -0.00324018}

I am having trouble integrating this, though:
F[v_] := NIntegrate[FxyzNew[x, y, z, r, 1/100], {r, -2, v}, {x, -Pi, Pi}, 
   {y, -Pi, Pi}, {z, -Pi, Pi}, Method -> {Automatic, "SymbolicProcessing" -> 0}];

Timing@F[1]

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 586.4907581840607` and 4625.0642323092015` for the integral and error estimates.

{1.90051, 586.491}

Looking at the integrand on a slice of parameter space
Plot[FxyzNew[0, y, 0.3, -1, 1/100], {y, -Pi, Pi}, PlotRange -> Full]

it seems like large cancellations need to occur in the integration.
However, I am very bad at numerical integration, so perhaps somebody has a good idea how to best handle a case like this.

How can I perform a successive NIntegrate process with two different variables?

One Answer

Add your own answers!

Ask a Question