How to numerical solve equation with 3x-integrals

Questions
, , ,
1

Good day everyone. I have a big problem with Mathematica. I cannot calculate equation with 3x-integrals. Equation will be needed solve numerically. Please, help me, I don’t know how to solve this.

I need to find Delta:

Delta = mu + 2*U*(sigma - rho)

where

rho = (1/2)*
  NIntegrate[(eps + Delta)/Ek - 1, {kx, -Pi, Pi}, {ky, -Pi, 
    Pi}, {kz, -Pi, Pi}]
sigma = -(1/2)*
  NIntegrate[Delta/Ek, {kx, -Pi, Pi}, {ky, -Pi, Pi}, {kz, -Pi, Pi}]

where

m = 0.02
J = 50
muB = 0.67
Dst = 7.1
U = 313
g = 2.06
mu = muB*g*Hext - Dst
eps = J*(3 - Cos[kx] - Cos[ky] - Cos[kz])
Ek = eps^(1/2)*(eps + Delta)^(1/2)

Here Hext should change from 6 to 12 with the step 0.5

Please, help me.

2

I get the impression that there is no solution of your problem.

And there are problems: may be your integrand has a singularity in the region of integration. I did not check this. I tried the following

m = 0.02
J = 50
muB = 0.67
Dst = 7.1
U = 313
g = 2.06
mu = muB*g*Hext - Dst
eps = J*(3 - Cos[kx] - Cos[ky] - Cos[kz])
Ek[Delta_] := eps^(1/2)*(eps + Delta)^(1/2)

(Note the error messages in what follows pointing to an unfriendly integrand. And it could be that the results are simply wrong)

trho = Table[
  (1/2)*NIntegrate[(eps + Delta)/Ek[Delta] - 1, {kx, -Pi, Pi}, {ky, -Pi, Pi}, {kz, -Pi, Pi}], {Delta, 0, 2, .2}]
tsigma = Table[
  -(1/2)*NIntegrate[Delta/Ek[Delta], {kx, -Pi, Pi}, {ky, -Pi, Pi}, {kz, -Pi, Pi}], {Delta, 0, 2, .2}]

And then

mu + 2*U*(tsigma - trho)[[1]] /. Hext -> 5.14418

Here you can see that Delta = 0 is a solution if Hext = 5.14418

3

Thank you for your answer. It will be helpful for me. I try to simplify the integral, previously. Thank you.

4

Ok. Your integrand has indeed a singularity at the origin. See below. Are you sure that your model / equations are correct?

m = 0.02
J = 50
muB = 0.67
Dst = 7.1
U = 313
g = 2.06
mu = muB*g*Hext - Dst
eps[kx_, ky_, kz_] := J*(3 - Cos[kx] - Cos[ky] - Cos[kz])
Ek[Delta_, kx_, ky_, kz_] := 
 eps[kx, ky, kz]^(1/2)*(eps[kx, ky, kz] + Delta)^(1/2)

eps[0, 0, 0]
Ek[del, 0, 0, 0]

So the expression containing Ek will tend to infinitiy whern Ek is evaulated at the origin

fac = .2;
Plot3D[((eps[kx, ky, 0] + 1)/Ek[2, kx, ky, 0] - 1) Boole[   kx^2 + ky^2 + 0 > .0004], {kx, -Pi fac, Pi fac}, {ky, -Pi fac, 
  Pi fac}, PlotRange -> All, PlotPoints -> 20]


Plot[((eps[kx, 0, 0] + 1)/Ek[2, kx, 0, 0] - 1), {kx, -Pi fac, Pi fac},
  PlotRange -> All, PlotStyle -> {Thick, Red}, AxesOrigin -> {0, 0}]
5

This I think re-inforces what Hans showed.

sol = Series[eq1, {kx, 0, 3}, {ky, 0, 3}, {kz, 0, 3}] // 
    Normal // ExpandAll;

For either eq1 or eq2 is a “typical” series with factors kz Delta in the denominator but sporadically powered so (i believe) neither can be made to eliminate the other in all terms of the series.

Limit[sol, kz → 0] ~= kx^(n-1) * ky^(n-1) * Sqrt[Delta] * Infinity

NIntegrate will approach inf when kz or Delta → 0. All assuming real values (Nintegrate does not handle complex values, and complex analysis of the problem would be much different).