# How to numerical solve equation with 3x-integrals

**URL:** <https://community.wolfram.com/t/how-to-numerical-solve-equation-with-3x-integrals/16790>\
**Category:** Questions\
**Tags:** algebra, calculus, equation-solving, numerical-computation\
**Created:** [September 28, 2020, 2:17pm UTC](https://community.wolfram.com/t/how-to-numerical-solve-equation-with-3x-integrals/16790 "2020-09-28T14:17:15Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![uzjavabiz](https://avatars.discourse-cdn.com/v4/letter/u/f6c823/32.png) [@uzjavabiz](https://community.wolfram.com/u/uzjavabiz)\
**Post date:** [September 28, 2020, 2:17pm UTC](https://community.wolfram.com/t/how-to-numerical-solve-equation-with-3x-integrals/16790/1 "2020-09-28T14:17:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![hdolhaine](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/hdolhaine/32/35924_2.png) [@hdolhaine](https://community.wolfram.com/u/hdolhaine)\
**Post date:** [September 28, 2020, 10:15pm UTC](https://community.wolfram.com/t/how-to-numerical-solve-equation-with-3x-integrals/16790/2 "2020-09-28T22:15:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![uzjavabiz](https://avatars.discourse-cdn.com/v4/letter/u/f6c823/32.png) [@uzjavabiz](https://community.wolfram.com/u/uzjavabiz)\
**Post date:** [September 29, 2020, 7:26pm UTC](https://community.wolfram.com/t/how-to-numerical-solve-equation-with-3x-integrals/16790/3 "2020-09-29T19:26:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![hdolhaine](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/hdolhaine/32/35924_2.png) [@hdolhaine](https://community.wolfram.com/u/hdolhaine)\
**Post date:** [September 29, 2020, 11:16pm UTC](https://community.wolfram.com/t/how-to-numerical-solve-equation-with-3x-integrals/16790/4 "2020-09-29T23:16:38Z")

</div>

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}]

```

---

<div class="post-metadata">

**Author:** ![anon\_91240867](https://avatars.discourse-cdn.com/v4/letter/a/8baadc/32.png) [@anon\_91240867](https://community.wolfram.com/u/anon_91240867)\
**Post date:** [September 30, 2020, 9:40am UTC](https://community.wolfram.com/t/how-to-numerical-solve-equation-with-3x-integrals/16790/5 "2020-09-30T09:40:37Z")

</div>

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).