Message Boards Message Boards

0
|
6366 Views
|
11 Replies
|
1 Total Likes
View groups...
Share
Share this post:
GROUPS:
Mathematica

How to propely define arbitrary constants?

Mitja Jan?i?, University
Posted 10 years ago

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!
POSTED BY: Mitja Jan?i?
11 Replies
Sort By:

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
POSTED BY: David Reiss

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
POSTED BY: Mitja Jan?i?

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.
POSTED BY: David Reiss

To Frank: I am calculating the Biot-Savart law (Equation 337 on http://farside.ph.utexas.edu/teaching/em/lectures/node39.html). 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!
POSTED BY: Mitja Jan?i?

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
POSTED BY: David Reiss

Ok, next question: Why are you only integrating the derivative of Koncna with respect to y, which is the second component of Elpolje ?
POSTED BY: Frank Kampas

It doesn't change anything. :D

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.
POSTED BY: Mitja Jan?i?

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?
POSTED BY: Frank Kampas

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.... -.-
POSTED BY: Mitja Jan?i?

To specify that your constants are Real, make use of the Assumptions option to Integrate.
POSTED BY: David Reiss 
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
POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract