Message Boards Message Boards

0
|
2001 Views
|
7 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Numerically solving an integral: no results

Posted 1 year ago

I am trying to numerically solve the following integral. However, I am not getting any result.

N[Integrate[(E^((-1/3000^2*(-3000 + x1)^2 - (-8000 + x2)^2/4000^2 - (-12000 + x3)^2/6000^2)/2) * (7.5*x1 + 6*(-x1 + x2) + 4.5*(-x2 + x3)))/(3000*4000*6000*2*Sqrt[2]*Pi^(3/2)),
  {x1,0,3000},{x2,x1,3000},{x3,x2,3000}]]
POSTED BY: Ma Wi
7 Replies
Posted 1 year ago

I think I have misunderstood your expression. Please edit your post and paste the InputForm of your expression with a blank line before and after and with four spaces at the beginning of the line. Or just attach your notebook file to your post.

Thank you.

POSTED BY: Bill Nelson
Posted 1 year ago

Thanks for your feedback. I edited it.

POSTED BY: Ma Wi
Posted 1 year ago

I am guessing that you wanted to change all

\(40)

to

(

and

\(41)

to

)

and

e

to

E

and

Integrate[...,{x3,x2,3000},{x2,x1,3000},{x1,0,3000}]

to

Integrate[...,{x1,0,3000},{x2,x1,3000},{x3,x2,3000}]

If all that is correct then your input becomes

Integrate[(E^((-1/9*(-3 + x1)^2 - (-8 + x2)^2/16 - (-12 + x3)^2/36)/2)*
  (7.5*x1 + 6*(-x1 + x2) + 4.5*(-x2 + x3)))/(144*Sqrt[2]*Pi^(3/2)),
  {x1,0,3000},{x2,x1,3000},{x3,x2,3000}]

Then this

Table[x1=RandomReal[{0,3000}];x2=RandomReal[{x1,3000}];x3=RandomReal[{x2,3000}];
  (E^((-1/9*(-3 + x1)^2 - (-8 + x2)^2/16 - (-12 + x3)^2/36)/2)*
  (7.5*x1 + 6*(-x1 + x2) + 4.5*(-x2 + x3)))/(144*Sqrt[2]*Pi^(3/2)),
  {1000}]

returns this

{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0....

seeming to show that your integrand is very small, almost zero, almost everywhere.

It appears that only for small values of x1,x2,x3, less than 50, is your integrand slightly larger than zero.That appears to be confusing the integration algorithm, where it looks at values of your integrand hundreds or thousands of times and finds the value is almost exactly zero every time.

If all of my changes and guesses are correct then this

NIntegrate[(E^((-1/9*(-3 + x1)^2 - (-8 + x2)^2/16 - (-12 + x3)^2/36)/2)*
  (7.5*x1 + 6*(-x1 + x2) + 4.5*(-x2 + x3)))/(144*Sqrt[2]*Pi^(3/2)),
  {x1,0,40},{x2,x1,40},{x3,x2,40}]

estimates your full integral is approximately 38.2249 but I do not know how much error is associated with that.

Please check every detail of this and try to make certain that you have found and corrected all my mistakes.

POSTED BY: Bill Nelson
Posted 1 year ago

Sorry for the confustion. It seems like the syntax of Wolfram Alpha and Mathematica is not completely compatible. All of your conversions are correct. Thanks for pointing out that the integrand is mostly zero. This was due to the fact that mean and standard deviation were specified in the wrong dimension. Multiplying all of them with a factor of 1000 yields the correct integral:

Integrate[(E^((-1/9000000*(-3000 + x1)^2 - (-8000 + x2)^2/16000000 - (-12000 + x3)^2/36000000)/2) * (7.5*x1 + 6*(-x1 + x2) + 4.5*(-x2 + x3)))/(3000*4000*6000*2*Sqrt[2]*Pi^(3/2)),
  {x1,0,3000},{x2,x1,3000},{x3,x2,3000}]

Unfortunately, I am still not getting any result.

POSTED BY: Ma Wi
Posted 1 year ago

I am sorry. I had no idea that you were trying to do this using WolframAlpha.

In WolframAlpha If I simplify your integral to

Integrate[(E^((-(-3000+x1)^2/9-(-8000+x2)^2/16-(-12000+x3)^2/36)/(2*10^6))*(1.5(x1+x2+3x3))),{x1,0,3000},{x2,x1,3000},{x3,x2,3000}]

WolframAlphacalculation

that barely completes within the time allowed for the free version of Wolfram Alpha and gives a result

3.90247*10^12.

For the final step I calculate

3.90247*10^12/(144*10^9*Sqrt[2]*Pi^(3/2))

and WolframAlpha gives a result of 3.44142 which equals the result from Wolfram Mathematica if I do the complete problem in one step.

Are you able to complete this calculation now?

Please check every step of this very carefully to find and fix any errors that I have made.

POSTED BY: Bill Nelson
Posted 1 year ago

Thanks a lot! With this simplification I was able to solve three out of my four integrals over the Gaussian distributions. For the final one, with the integrant being 1, I do not receive any result with the exact upper bounds being infinity:

Integrate[E^((-(-3000+x1)^2/9-(-8000+x2)^2/16-(-12000+x3)^2/36)/(2*10^6)),{x1,3000,Infty},{x2,x1,Infty},{x3,x2,Infty}]

So I tried setting the upper bounds to sufficiently high values such as

Integrate[E^((-(-3000+x1)^2/9-(-8000+x2)^2/16-(-12000+x3)^2/36)/(2*10^6)),{x1,3000,200000},{x2,x1,500000},{x3,x2,500000}]

and I obtain 9.25339x10^9. However, it confuses me that if I set the upper bounds of x2 and x3 from 5x10^5 to 6x10^5, I am getting a quite different result, namely 1.79312x10^8.

Integrate[E^((-(-3000+x1)^2/9-(-8000+x2)^2/16-(-12000+x3)^2/36)/(2*10^6)),{x1,3000,200000},{x2,x1,600000},{x3,x2,600000}]

Moreover, I thought this might be due to numerical instabilities because of the large values, so I rescaled all xi by 1/1000:

Integrate[10^9*E^((-(-3+x1)^2/9-(-8+x2)^2/16-(-12+x3)^2/36)/(2)),{x1,3,200},{x2,x1,500},{x3,x2,500}]

and I obtain the same result as before the scaling, 9.25339x10^9. And the same for the upper bounds of x2 and x3 set to 600: 1.79312x10^8. How can the result get smaller when extending the range of the integration over a (very small) positive function?

POSTED BY: Ma Wi

The integral from the question/first post takes only a few milliseconds with NIntegrate[]:

NIntegrate[(E^((-1/3000^2*(-3000 + x1)^2 - (-8000 + x2)^2/
          4000^2 - (-12000 + x3)^2/6000^2)/2)*(7.5*x1 + 
      6*(-x1 + x2) + 4.5*(-x2 + x3)))/(3000*4000*6000*2*Sqrt[2]*
    Pi^(3/2)), {x1, 0, 3000}, {x2, x1, 3000}, {x3, x2, 3000}]

(*  3.44142  *)
POSTED BY: Michael Rogers
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