I have one really horrible integral to calculate and Mathematica doesn’t want to do it.
I suspect the problems are my constants (because if I leave them out, it does everything perfect but I need them). So how do I properly define constants?
They are real numbers, not complex.
For example in the integral is also expression: |r-r0| , where r and r0 are both vectors but r0 is a constant vector.
Thank you for help!
In[1]:= r0 = {1, 2, 3};
In[5]:= NIntegrate[
Norm[{x, y, z} - r0], {x, -5, 5}, {y, -5, 5}, {z, -5, 5}]
Out[5]= 5870.68
To specify that your constants are Real, make use of the Assumptions option to Integrate.
To Frank: I didn’t like your idea, because in the end I would like to have the result: B(x0,y0,z0) if my x0,y0,z0 are the constants I am talking about, but I tried it anyway. (see below).
To David: Yours doesn’t work either.
I still get a complex result. Here is the code:
In[2]:= Elpolje = Simplify[-Grad[Koncna[x, y], {x, y}]]
Out[2]= {(
4 a A (a^2 - x^2 + y^2))/((a^2 - 2 a x + x^2 + y^2) (a^2 + 2 a x +
x^2 + y^2)), -((
8 a A x y)/((a^2 - 2 a x + x^2 + y^2) (a^2 + 2 a x + x^2 + y^2)))}
The integral
In[8]:= a = 20
y0 = 1
x0 = 1
z0 = 1
Out[8]= 20
Out[9]= 1
Out[10]= 1
Out[11]= 1
In[12]:= Integrate[(Elpolje[[2]]*
z0)/((aA) ((x0 - x)^2 + (y0 - y)^2 + (z0 - z)^2)^(3/2)), x]
Out[12]= -(1/aA)
160 A y ((-482379 + 157585 x + 326374 y + 1628 x y - 168785 y^2 -
20 x y^2 + 4796 y^3 + 8 x y^3 - 1600 y^4 + 326374 z +
1628 x z - 6408 y z - 24 x y z + 4812 y^2 z + 8 x y^2 z -
8 y^3 z - 166391 z^2 - 826 x z^2 + 3228 y z^2 + 12 x y z^2 -
2410 y^2 z^2 - 4 x y^2 z^2 + 4 y^3 z^2 + 3212 z^3 + 12 x z^3 -
24 y z^3 + 4 y^2 z^3 - 813 z^4 - 3 x z^4 + 6 y z^4 - y^2 z^4 +
6 z^5 - z^6)/(Sqrt[
3 - 2 x + x^2 - 2 y + y^2 - 2 z +
z^2] (51719068962 - 52755965394 y + 27935849089 y^2 -
1568132376 y^3 + 535095200 y^4 - 10253120 y^5 +
2560064 y^6 - 52755965394 z + 2089332928 y z -
1583724024 y^2 z + 28405040 y^3 z - 10330368 y^4 z +
25728 y^5 z + 27422649161 z^2 - 1068053936 y z^2 +
806167828 y^2 z^2 - 14318520 y^3 z^2 + 5178112 y^4 z^2 -
12864 y^5 z^2 - 1052462288 z^3 + 23490768 y z^3 -
14344696 y^2 z^3 + 116064 y^3 z^3 - 12928 y^4 z^3 +
272886176 z^4 - 6001972 y z^4 + 3634790 y^2 z^4 -
29096 y^3 z^4 + 3232 y^4 z^4 - 5898548 z^5 + 77792 y z^5 -
29192 y^2 z^5 + 48 y^3 z^5 + 1013294 z^6 - 13152 y z^6 +
4884 y^2 z^6 - 8 y^3 z^6 - 12992 z^7 + 80 y z^7 -
8 y^2 z^7 + 1654 z^8 - 10 y z^8 + y^2 z^8 - 10 z^9 + z^10)) +
Log[-((320 y (-363 I + (38 + 2 I) y + 2 I z - I z^2) Sqrt[
3 - 2 x + x^2 - 2 y + y^2 - 2 z + z^2])/(-20 + x - I y)) + (
320 (6171 y -
6897 x y + (692 + 1009 I) y^2 + (38 - 1085 I) x y^2 - (405 -
74 I) y^3 + (38 + 2 I) x y^3 + (2 - 38 I) y^4 + 692 y z +
38 x y z - (8 - 74 I) y^2 z + 2 I x y^2 z + 2 y^3 z -
350 y z^2 - 19 x y z^2 + (4 - 37 I) y^2 z^2 - I x y^2 z^2 -
y^3 z^2 + 4 y z^3 - y z^4))/((-20 I + I x + y) Sqrt[
363 - (2 - 38 I) y - 2 z + z^2])]/(
160 y (-363 I + (38 + 2 I) y + 2 I z - I z^2) Sqrt[
363 - (2 - 38 I) y - 2 z + z^2]) + (
I Log[(
320 I y (-443 + (2 + 42 I) y + 2 z - z^2) Sqrt[
3 - 2 x + x^2 - 2 y + y^2 - 2 z + z^2])/(20 + x - I y) - (
320 (-10189 y +
9303 x y + (932 + 1409 I) y^2 - (42 +
1325 I) x y^2 - (405 + 86 I) y^3 - (42 -
2 I) x y^3 + (2 + 42 I) y^4 + 932 y z -
42 x y z - (8 + 86 I) y^2 z + 2 I x y^2 z + 2 y^3 z -
470 y z^2 + 21 x y z^2 + (4 + 43 I) y^2 z^2 - I x y^2 z^2 -
y^3 z^2 + 4 y z^3 - y z^4))/((20 I + I x + y) Sqrt[
443 - (2 + 42 I) y - 2 z + z^2])])/(
160 y (-443 + (2 + 42 I) y + 2 z - z^2) Sqrt[
443 - (2 + 42 I) y - 2 z + z^2]) + (
I Log[(320 I y (-363 + (2 + 38 I) y + 2 z - z^2) Sqrt[
3 - 2 x + x^2 - 2 y + y^2 - 2 z + z^2])/(-20 + x + I y) + (
320 (6171 y -
6897 x y + (692 - 1009 I) y^2 + (38 +
1085 I) x y^2 - (405 + 74 I) y^3 + (38 -
2 I) x y^3 + (2 + 38 I) y^4 + 692 y z +
38 x y z - (8 + 74 I) y^2 z - 2 I x y^2 z + 2 y^3 z -
350 y z^2 - 19 x y z^2 + (4 + 37 I) y^2 z^2 + I x y^2 z^2 -
y^3 z^2 + 4 y z^3 - y z^4))/((20 I - I x + y) Sqrt[
363 - (2 + 38 I) y - 2 z + z^2])])/(
160 y (-363 + (2 + 38 I) y + 2 z - z^2) Sqrt[
363 - (2 + 38 I) y - 2 z + z^2]) +
Log[-((320 y (-443 I + (42 + 2 I) y + 2 I z - I z^2) Sqrt[
3 - 2 x + x^2 - 2 y + y^2 - 2 z + z^2])/(20 + x + I y)) + (
320 (10189 y -
9303 x y - (932 - 1409 I) y^2 + (42 - 1325 I) x y^2 + (405 -
86 I) y^3 + (42 + 2 I) x y^3 - (2 - 42 I) y^4 - 932 y z +
42 x y z + (8 - 86 I) y^2 z + 2 I x y^2 z - 2 y^3 z +
470 y z^2 - 21 x y z^2 - (4 - 43 I) y^2 z^2 - I x y^2 z^2 +
y^3 z^2 - 4 y z^3 + y z^4))/((-20 I - I x + y) Sqrt[
443 - (2 - 42 I) y - 2 z + z^2])]/(
160 y (-443 I + (42 + 2 I) y + 2 I z - I z^2) Sqrt[
443 - (2 - 42 I) y - 2 z + z^2]))
I ignored the whole result as soon as I say those imaginary parts in it… -.-
I was only trying to show how to handle the constant vector. Why are you only integrating over x in your last post?
Don’t you want to integrate over y and z also?
It doesn’t change anything. 
If the integration over x is wrong than everything is wrong. That’s one reason, and the other reason is that integrating this result above also over y is (I guess) to complex for mathematica and I can wait for 5 minutes but still nothing happens so I am forced to abort the process.
AND I am highly suspecting that those complex numbers are the main reason, why Mathematica can’t do the second integration over y.
Ok, next question: Why are you only integrating the derivative of Koncna with respect to y, which is the second component of Elpolje ?
In the course of doing the integral Mathematica is keeping track of the various potential branch cuts. The end result will nonetheless be real. You should check things numerically to verify this.
For example integrate the following:
temp = Integrate[Sqrt[ x]/(1 + (x - 1)^2), x]
this gives
-((2 ArcTan[Sqrt[x]/Sqrt[-1 - I]])/(-1 - I)^(3/2)) -
(2 ArcTan[Sqrt[x]/Sqrt[-1 + I]])/(-1 + I)^(3/2)
For a particular value of x, say the following,
N[temp] /. x -> 2
it gives
1.4994 + 0. I
In fact, in this particular case
-((2 ArcTan[Sqrt[x]/Sqrt[-1 - I]])/(-1 - I)^(3/2)) -
(2 ArcTan[Sqrt[x]/Sqrt[-1 + I]])/(-1 + I)^(3/2)
is the sum of an expression and its complex conjugate.
I suggest doing some numerical experiments with your first integration (using particular values of a, A, x0, etc… to confirm this
To Frank:
I am calculating the Biot-Savart law (Equation 337 on The Biot-Savart law). Originaly there should also be the third component in the nominator, but the third component of Elpolje is obviously 0.
To David:
Yeah that’s why I know that something is wrong. Everything is real, so should be the result!
But the result is real.
Here is your integral (the one only with respect to x that you show above and where I fix the aA to a A – i.e., add a space):
temp = Integrate[((-((
8 a A x y)/((a^2 - 2 a x + x^2 + y^2) (a^2 + 2 a x + x^2 +
y^2))))*
z0)/((a A) ((x0 - x)^2 + (y0 - y)^2 + (z0 - z)^2)^(3/2)), x];
now test it with random values for the unspecified parameters:
In[95]:= Table[
Im[temp /. {x0 -> Random[], y0 -> Random[], z0 -> Random[],
A -> Random[], a -> Random[], y -> Random[], z -> Random[],
x -> Random[]} // N], {1000}] // Union
Out[95]= {0.}
The imaginary part is always 0.
Yes, true. My bad!
Hmm, well integration from - infinity to infinity over x and y simply doesn’t work. So I switched to numerical integration over exactly defined area. But I get Impossible errors that I simply DO NOT UNDERSTAND
In[3]:= Elpolje
Out[3]= {(
4 a A (a^2 - x^2 + y^2))/((a^2 - 2 a x + x^2 + y^2) (a^2 + 2 a x +
x^2 + y^2)), -((
8 a A x y)/((a^2 - 2 a x + x^2 + y^2) (a^2 + 2 a x + x^2 + y^2)))}
In[60]:= x0 = 10
y0 = 10
z0 = 10
a = 20
Out[60]= 10
Out[61]= 10
Out[62]= 10
Out[63]= 20
In[64]:= NIntegrate[
NIntegrate[(Elpolje[[2]]*z0)/(a*
A*((x0 - x)^2 + (y0 - y)^2 + z0^2)^(3/2)), {x, -1000,
1000}], {y, -1000, 1000}]
During evaluation of In[64]:= NIntegrate::inumr: The integrand -((80 x y)/((980200+(10-x)^2)^(3/2) (400-40 x+x^2+y^2) (400+40 x+x^2+y^2))) has evaluated to non-numerical values for all sampling points in the region with boundaries {{-1000,0.}}. >>
During evaluation of In[64]:= NIntegrate::inumr: The integrand -((80 x y)/((980200+(10-x)^2)^(3/2) (400-40 x+x^2+y^2) (400+40 x+x^2+y^2))) has evaluated to non-numerical values for all sampling points in the region with boundaries {{-1000,0.}}. >>
During evaluation of In[64]:= NIntegrate::inumr: The integrand -((80 x y)/((980200+(10-x)^2)^(3/2) (400-40 x+x^2+y^2) (400+40 x+x^2+y^2))) has evaluated to non-numerical values for all sampling points in the region with boundaries {{-1000,0.}}. >>
During evaluation of In[64]:= General::stop: Further output of NIntegrate::inumr will be suppressed during this calculation. >>
During evaluation of In[64]:= 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. >>
During evaluation of In[64]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in y near {y} = {3.87517}. NIntegrate obtained -0.00901131 and 0.00018888873826197088` for the integral and error estimates. >>
Out[64]= -0.00901131
This is happening because you have an NIntegrate within another NIntegrate. The inner one has the variable y as a nonnumerical parameter.
You want to use a single NIntegrate with both integration variables – x and y – as iterators. This is the form you want:
NIntegrate[f, {x, xmin, xmax}, {y, ymin, ymax}]
or
NIntegrate[f, {y, ymin, ymax}, {x, xmin, xmax}]
Also you have not given A a value.
Here is your expression with these things fixed.:
In[6]:=
x0 = 10;
y0 = 10;
z0 = 10;
a = 20;
A = 1;
In[12]:=
NIntegrate[((-((8 a A x y)/((a^2 - 2 a x + x^2 + y^2) (a^2 + 2 a x +
x^2 + y^2))))*z0)/(a*
A*((x0 - x)^2 + (y0 - y)^2 + z0^2)^(3/2)), {x, -1000,
1000}, {y, -1000, 1000}]
Out[12]= -0.00972789