Message Boards Message Boards

Integration involving vectors?

Posted 7 years ago

I am trying to integrate the following expression - where ${\bf p3}$, ${\bf q3}$, ${\bf p}$, ${\bf q}$ are vectors. The angle between each pair of vectors e.g. ${\bf q3}$, ${\bf p3}$ can't be ignored. $\mu$ is a scalar.

$\int d{\bf p3} \int d{\bf q3} \int d{\bf q} \frac{1}{\sqrt{\mu ^2+\text{q3}^2} \sqrt{\mu ^2+({\bf p}-{\bf p3}-{\bf q3})^2} \sqrt{\mu ^2+(-{\bf p3}+{\bf q}-{\bf q3})^2} \left(\sqrt{\mu ^2+(-{\bf p3}+{\bf q}-{\bf q3})^2}+\sqrt{\mu ^2+{\bf q3}^2}\right)}$.

Tried the following:

q3v = Table[Subscript[q3, i], {i, 3}]; 
p3v = Table[Subscript[p3, i], {i, 3}];
qv = Table[Subscript[q, i], {i, 3}];
 pv = Table[Subscript[p, i], {i, 3}]; 

and then

a := 1/Sqrt[q3v^2 + \[Mu]^2]; 
b := 1/(Sqrt[q3v^2 + \[Mu]^2] + Sqrt[(-p3v + qv - q3v)^2 + \[Mu]^2]); 
c := 1/(Sqrt[(pv - p3v - q3v)^2 + \[Mu]^2] Sqrt[(-p3v + qv - 
        q3v)^2 + \[Mu]^2])

The output of a.b.c gives

{1/(Sqrt[\[Mu]^2 + (Subscript[p, 1] - Subscript[p3, 1] - Subscript[q3, 1])^2]*Sqrt[\[Mu]^2 + (-Subscript[p3, 1] + Subscript[q, 1] - Subscript[q3, 1])^2]*Sqrt[\[Mu]^2 + Subscript[q3, 1]^2]*
    (Sqrt[\[Mu]^2 + (-Subscript[p3, 1] + Subscript[q, 1] - Subscript[q3, 1])^2] + Sqrt[\[Mu]^2 + Subscript[q3, 1]^2])), 
  1/(Sqrt[\[Mu]^2 + (Subscript[p, 2] - Subscript[p3, 2] - Subscript[q3, 2])^2]*Sqrt[\[Mu]^2 + (-Subscript[p3, 2] + Subscript[q, 2] - Subscript[q3, 2])^2]*Sqrt[\[Mu]^2 + Subscript[q3, 2]^2]*
    (Sqrt[\[Mu]^2 + (-Subscript[p3, 2] + Subscript[q, 2] - Subscript[q3, 2])^2] + Sqrt[\[Mu]^2 + Subscript[q3, 2]^2])), 
  1/(Sqrt[\[Mu]^2 + (Subscript[p, 3] - Subscript[p3, 3] - Subscript[q3, 3])^2]*Sqrt[\[Mu]^2 + (-Subscript[p3, 3] + Subscript[q, 3] - Subscript[q3, 3])^2]*Sqrt[\[Mu]^2 + Subscript[q3, 3]^2]*
    (Sqrt[\[Mu]^2 + (-Subscript[p3, 3] + Subscript[q, 3] - Subscript[q3, 3])^2] + Sqrt[\[Mu]^2 + Subscript[q3, 3]^2]))}

Now, let's try to evaluate the first of the triplet in the output of a.b.c

Integrate[1/(Sqrt[\[Mu]^2 + (Subscript[p, 1] - Subscript[p3, 1] - Subscript[q3, 1])^2]*Sqrt[\[Mu]^2 + (-Subscript[p3, 1] + Subscript[q, 1] - Subscript[q3, 1])^2]*Sqrt[\[Mu]^2 + Subscript[q3, 1]^2]*
    (Sqrt[\[Mu]^2 + (-Subscript[p3, 1] + Subscript[q, 1] - Subscript[q3, 1])^2] + Sqrt[\[Mu]^2 + Subscript[q3, 1]^2])), {q3, -1, 1}]

which Mathematica 11 will not evaluate.

Is this approach correct?

POSTED BY: Arny Toynbee
13 Replies
Posted 7 years ago

Thank you, Hans.

POSTED BY: Arny Toynbee

jjm = jj /. m -> 2. implies that m is being set to 2. prior to performing the integration.

Yes. To do the numerical integration m must have a definite (numerical) value.

is it over all of p3v, q3v, pv, qv?

Yes. All p3v, q3v, pv, qv? Look at intvariables2 to see the boundaries. See attached notebook.

Attachments:
POSTED BY: Hans Dolhaine

At last I noticed that ther is a typo in the definition of the integrand. It should read (at least I think so)

a := 1/Sqrt[m^2 + q3v.q3v]
b := 1/Sqrt[m^2 + (p3v + q3v - pv).(p3v + q3v - pv)]
c := 1/Sqrt[m^2 + (p3v + q3v - qv).(p3v + q3v - qv)]
d := 1/(Sqrt[m^2 + (p3v + q3v - qv).(p3v + q3v - qv)] + Sqrt[m^2 + q3v.q3v])

With this

In[9]:= jj = a b c d // FullSimplify

Out[9]= 1/(Sqrt[
   m^2 + q3[1]^2 + q3[2]^2 +  q3[3]^2] \[Sqrt](m^2 + (-p[1] + p3[1] + q3[1])^2 + 
  (-p[2] + p3[2] + q3[2])^2 + (-p[3] + p3[3] + q3[3])^2) Sqrt[m^2 + (p3[1] - q[1] + q3[1])^2 + (p3[2] - q[2] + q3[2])^2+
   (p3[3] - q[3] + q3[3])^2] (Sqrt[ m^2 + q3[1]^2 + q3[2]^2 + q3[3]^2] + \[Sqrt](m^2 + (p3[1] - q[1] + q3[1])^2 +
   (p3[2] - q[2] + q3[2])^2 + (p3[3] - q[3] + q3[3])^2)))

and the denominator should never vanish.

Then

jjm = jj /. m -> 2.

and

In[15]:= NIntegrate[jjm, Evaluate[Sequence @@ intvariables2], 
 MaxPoints -> 100]

During evaluation of In[15]:= NIntegrate::maxp: The integral failed to converge after 13227 integrand evaluations. NIntegrate obtained 60.37299593570131` and 8.705681178065333` for the integral and error estimates. >>

Out[15]= 60.373

The Integration is reasonably fast. Without the MaxPoint-Option it takes quite a time, so I aborted it. But it could be an idea to let it run and compare the result with that with the MaxPoint-Option.

POSTED BY: Hans Dolhaine
Posted 7 years ago

Thank you. Sorry for the following really dumb question - jjm = jj /. m -> 2. implies that m is being set to 2. prior to performing the integration. What is the @@ intvariables2 in this case, is it over all of p3v, q3v, pv, qv?

POSTED BY: Arny Toynbee

First of all some remarks. I am pretty sure that your code gives a wrong integrand.

Your formula involves the square of vectors. But these are given by dot products. If you do it as you do you get a list of squares, and my ( I write m ) added to it gives another list or vector with m added to each component, because Plus is distributed over lists:

In[1]:= v = {a, b, c};
v^2
m + v^2

Out[2]= {a^2, b^2, c^2}

Out[3]= {a^2 + m, b^2 + m, c^2 + m}

But actually you can't add a scalar to a vector, these types don't fit together. Mathematica does it because the functions used are listable.

What about the squareroot? Sqrt is also distributed and so you get a vector (or list) of squareroots:

In[4]:= Sqrt[m + v^2]

Out[4]= {Sqrt[a^2 + m], Sqrt[b^2 + m], Sqrt[c^2 + m]}

And you use this in an expression as denominator with numerator 1. This gives a vector of fractions:

In[5]:= 1/Sqrt[m + v^2]

Out[5]= {1/Sqrt[a^2 + m], 1/Sqrt[b^2 + m], 1/Sqrt[c^2 + m]}

I am pretty sure the real interpretation of your formula is this

In[6]:= 1/Sqrt[m^2 + v.v]

Out[6]= 1/Sqrt[a^2 + b^2 + c^2 + m^2]

Having said this I redefine your integrand (note that I omit the underscripts. It was often remarked in this forum that variables with underscripts may behave sometimes in a faulty way)

    q3v = Table[q3[i], {i, 3}]
    p3v = Table[p3[i], {i, 3}]
    qv = Table[q[i], {i, 3}]
    pv = Table[p[i], {i, 3}]

    Out[10]= {q3[1], q3[2], q3[3]}

    Out[11]= {p3[1], p3[2], p3[3]}

    Out[12]= {q[1], q[2], q[3]}

Out[13]= {p[1], p[2], p[3]}


a := 1/Sqrt[m^2 + q3v.q3v]

b := 1/Sqrt[m^2 + (pv - p3v - q3v).(pv - p3v - q3v)]

c := 1/Sqrt[m^2 + (p3v - qv + q3v).(p3v + q3v + qv)]

d := 1/(Sqrt[m^2 + (p3v + q3v - qv).(p3v + q3v - qv)] +     Sqrt[m^2 + q3v.q3v])

It seems always to be ( p3v + q3v plus or minus something ) if you take minus-signs out of the brackets.

With this you have as integrand

In[19]:= jj = a b c d // FullSimplify

Out[19]= 1/(Sqrt[
   m^2 + q3[1]^2 + q3[2]^2 + 
    q3[3]^2] \[Sqrt](m^2 - q[1]^2 - q[2]^2 - 
      q[3]^2 + (p3[1] + q3[1])^2 + (p3[2] + q3[2])^2 + (p3[3] + 
        q3[3])^2) \[Sqrt](m^2 + (-p[1] + p3[1] + q3[1])^2 + (-p[2] + 
        p3[2] + q3[2])^2 + (-p[3] + p3[3] + q3[3])^2) (Sqrt[
     m^2 + q3[1]^2 + q3[2]^2 + 
      q3[3]^2] + \[Sqrt](m^2 + (p3[1] - q[1] + q3[1])^2 + (p3[2] - 
          q[2] + q3[2])^2 + (p3[3] - q[3] + q3[3])^2)))

ok. Now for the integration variables. I write

In[49]:= intvariables = Flatten[Join[{q3v}, {p3v}, {qv}, {pv}]]

Out[49]= {q3[1], q3[2], q3[3], p3[1], p3[2], p3[3], q[1], q[2], q[3],  p[1], p[2], p[3]}

and then

Integrate[jj, intvariables]

I had run this for about two hours without answer. It seems that this twelvefold integration is too complicated to be done symbolically as Frank remarked above.

However, you can do it numerically. Define the integrand by specifying my

jjm = jj /. m -> 2.

and fix some Integration limits (I do it here for all variables between 1 and 2. You may of course choose others for each variable

In[20]:= intvariables2 = {#, 1, 2} & /@ intvariables

Out[20]= {{q3[1], 1, 2}, {q3[2], 1, 2}, {q3[3], 1, 2}, {p3[1], 1,  2}, {p3[2], 1, 2}, {p3[3], 1, 2}, {q[1], 1, 2}, {q[2], 1, 2}, 
{q[3],1, 2}, {p[1], 1, 2}, {p[2], 1, 2}, {p[3], 1, 2}}

Then

In[22]:= NIntegrate[jjm, Evaluate[Sequence @@ intvariables2]]

Out[22]= 0.00292856
POSTED BY: Hans Dolhaine
Posted 7 years ago

Thank you for the corrections, and the caution on using subscripts. I have suffered for a long time with subscripts, without realizing this potential bug.

NIntegrate between -1, 1 produces the following output:\

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

and after 10+ minutes\

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 441.218 -211.434 I and 1135.6427949097094` for the integral and error estimates.\

441.218293633602 - 211.43356703452153*I

Is there a way in Mathematica to analyze the integrand behaviour in a way that can catch any ill-behaved regions?

Thanks

POSTED BY: Arny Toynbee

I don't know.

I got ( integrating for all variables from -1 to 1 )

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 70.76236375923493` and 0.010212021932754334` for the integral and error estimates. >>

With the above code I did

(# -> 0 & /@ intvariables)
jjm /. Drop[(# -> 0 & /@ intvariables), {1}]
jjm /. Drop[(# -> 0 & /@ intvariables), {5}]
jjm /. Drop[(# -> 0 & /@ intvariables), {12}]

and

Plot[jjm /. Drop[(# -> 0 & /@ intvariables), {1}], {q3[1], -1, 1}]
Plot[jjm /. Drop[(# -> 0 & /@ intvariables), {5}], {p3[2], -1, 1}]
Plot[jjm /. Drop[(# -> 0 & /@ intvariables), {12}], {p[3], -1, 1}]

to plot the integrand for some examples with all but one variable = zero. Seems to show a reasonable behavior.

then I tried

NIntegrate[jj /. m -> N[2., 30], Evaluate[Sequence @@ intvariables2], WorkingPrecision -> 30]

but this is still running (I didn't notice but I think for more than an hour)

POSTED BY: Hans Dolhaine

At last I noticed that ther is a typo in the definition of the integrand. It should read (at least I think so)

a := 1/Sqrt[m^2 + q3v.q3v]
b := 1/Sqrt[m^2 + (p3v + q3v - pv).(p3v + q3v - pv)]
c := 1/Sqrt[m^2 + (p3v + q3v - qv).(p3v + q3v - qv)]
d := 1/(Sqrt[m^2 + (p3v + q3v - qv).(p3v + q3v - qv)] + Sqrt[m^2 + q3v.q3v])

With this

In[9]:= jj = a b c d // FullSimplify

Out[9]= 1/(Sqrt[
   m^2 + q3[1]^2 + q3[2]^2 +  q3[3]^2] \[Sqrt](m^2 + (-p[1] + p3[1] + q3[1])^2 + 
  (-p[2] + p3[2] + q3[2])^2 + (-p[3] + p3[3] + q3[3])^2) Sqrt[m^2 + (p3[1] - q[1] + q3[1])^2 + (p3[2] - q[2] + q3[2])^2+
   (p3[3] - q[3] + q3[3])^2] (Sqrt[ m^2 + q3[1]^2 + q3[2]^2 + q3[3]^2] + \[Sqrt](m^2 + (p3[1] - q[1] + q3[1])^2 +
   (p3[2] - q[2] + q3[2])^2 + (p3[3] - q[3] + q3[3])^2)))

and the denominator should never vanish.

Then

jjm = jj /. m -> 2.

and

In[15]:= NIntegrate[jjm, Evaluate[Sequence @@ intvariables2], 
 MaxPoints -> 100]

During evaluation of In[15]:= NIntegrate::maxp: The integral failed to converge after 13227 integrand evaluations. NIntegrate obtained 60.37299593570131` and 8.705681178065333` for the integral and error estimates. >>

Out[15]= 60.373

The Integration is reasonably fast. Without the MaxPoint-Option it takes quite a time, so I aborted it. But it could be an idea to let it run and compare the result with that with the MaxPoint-Option.

POSTED BY: Hans Dolhaine

suggest you try a simpler version of the problem which will be easier to troubleshoot

POSTED BY: Frank Kampas
Posted 7 years ago

NIntegrate doesn't seem to work either.

POSTED BY: Arny Toynbee

You could do numerical integration for various values of lambda.

POSTED BY: Frank Kampas

Assuming your code is correct, it's probably too complicated to be integrated symbolically.

POSTED BY: Frank Kampas
Posted 7 years ago

Maple 17 appears to give an output when integrated wrt q3 for example, but the answer is very (several pages print) long. Is there a way you would suggest to integrate numerically, say between 0 and \lambda ? Thanks.

POSTED BY: Arny Toynbee
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract