Another Bessel function integral

zetafunc · 2017-02-21 10:06:27

The integral I would like to compute numerically is:

which, when written in terms of all of its vector components, looks like this:

where denotes the Bessel function of the first kind.

A sample Mathematica code would be

NIntegrate[(NIntegrate[(BesselJ[1, Sqrt[x^2 + y^2]]*
       BesselJ[1, Sqrt[(b - x)^2 + (c - y)^2]])/(Sqrt[x^2 + y^2]*
       Sqrt[(b - x)^2 + (c - y)^2]), {x, 0, 10}, {y, 0, 10}])^2, {b, 
  0, 10}, {c, 0, 10}]

but this returns something which does not make sense (it is a bunch of text involving LevinRuleDump and other similar things).

Some observations which can be made about this integral:

-It is bounded and continuous, and in particular is integrable in a neighbourhood of the origin. Therefore, any divergence can only come from integrating away from the origin.
-The 1-dimensional version of this integral converges numerically at infinity (that is, take away any of the integrals involving y and c, and set y = c = 0 in the integral -- this converges numerically, at least when integrated over [-10^100, 10^100]).
-The integral only converges conditionally, not absolutely (if it converges at all, which in principle it should).

Can anyone help find an estimate for this integral?

bobbym · 2017-02-21 11:02:58

Hi;

There seems to be some error in there somewhere. I am looking at a reformulation.

zetafunc · 2017-02-21 11:06:01

OK, let me know what you can find.

bobbym · 2017-02-21 11:14:52

It looks intractable. Usually, one has to resort to Monte Carlo methods for multidimensional problems. But I will try...

zetafunc · 2017-02-21 11:22:45

Would those methods be applicable here?

bobbym · 2017-02-21 12:34:58

I am not an expert on any of those and it is one of the reasons I stay away from 3D+.

To get useful ideas about any integration problem one has to look at the graph. Then it will be clear how to proceed. Take my integral as an example:

It is extremely tough but because I can graph the integrand I can gain insight into possble methods to deal with it. Yours is a 4 dimensional integral. I have never seen a 4D graph in our physical reality. The standard methods they teach in school are very weak. Even if we get a MCM working it will be very computationally intensive and will only yield 2 decimal places or so. This may not be enough to determine convergence.

Let me see what can be done and I will post my conclusions.

zetafunc · 2017-02-21 12:44:36

I would be very interested to hear your thoughts.

If the square causes a problem, there is a 6-integral representation also. I am currently trying this out with Mathematica. Rewriting it like this could help:

NIntegrate[(BesselJ[1, Sqrt[x^2 + y^2]]*
    BesselJ[1, Sqrt[z^2 + w^2]] BesselJ[1, 
     Sqrt[(b - x)^2 + (c - y)^2]]*
    BesselJ[1, Sqrt[(b + z)^2 + (c + w)^2]])/(Sqrt[x^2 + y^2]*
    Sqrt[z^2 + w^2]*Sqrt[(b - x)^2 + (c - y)^2]*
    Sqrt[(b + z)^2 + (c + w)^2]), {x, 0, 5}, {y, 0, 5}, {b, 0, 
  5}, {c, 0, 5}, {z, 0, 5}, {w, 0, 5}]

I am currently running that with those limits.

Last edited by zetafunc (2017-02-21 12:52:47)

zetafunc · 2017-02-21 13:05:16

I seem to be getting the result 55.1609 and 598.9347380610473 for the integral and error estimates, integrating from 0 to 10^6. The computation took 125.55 seconds. I used:

NIntegrate[(BesselJ[1, Sqrt[x^2 + y^2]]*
    BesselJ[1, Sqrt[z^2 + w^2]] BesselJ[1, 
     Sqrt[(b - x)^2 + (c - y)^2]]*
    BesselJ[1, Sqrt[(b + z)^2 + (c + w)^2]])/(Sqrt[x^2 + y^2]*
    Sqrt[z^2 + w^2]*Sqrt[(b - x)^2 + (c - y)^2]*
    Sqrt[(b + z)^2 + (c + w)^2]), {x, 0, 1000000}, {y, 0, 
  1000000}, {b, 0, 1000000}, {c, 0, 1000000}, {z, 0, 
  1000000}, {w, 0, 1000000}]

I've managed to do the computation up to 10^21. The trouble is, the error is huge (about 20000). However, the error does get smaller if we integrate over a larger and larger region. For instance, integrating up to 10^30 gives -506 with an error of 5422. Up t0 10^100, we get 121.589 and 392588.3228022117 for the integral and error estimates.

Last edited by zetafunc (2017-02-21 13:24:40)

bobbym · 2017-02-21 15:43:15

Hi;

I am getting an error almost as big as the answer with other messages about slow convergence. There is no reliability in that answer at all.

We have talked quite a bit about lots of stuff. I do not care about your physical appearance, your gender, your religion, your social status or how much money you have. I only know you from your writings, so in the truest sense I consider you a friend. So as your friend I would like to offer some advice:

Here is a quote from someone over at Mathematica.SE for an integral that is not as complex as yours.

Felix wrote:

You have lots of stuff in the integral that does not depend on x. You should first manually reduce the complexity of the integrand and the number of variables to the minimum possible. And then you need to use assumptions for all remaining parameters since the integral most likely does not converge for arbitrary complex values of the parameters. Even then it might not work, but in the way you write it even much simpler integrals would not give a result.

It is the job of the human to prepare the problem for M to solve. That is our job now. Mathematica can not do the impossible. I have watched you suffering over these type of multiple integrals for months with little progress being made. I know you can make progress when the problem allows it. This type of problem is not amenable to progress either by you or M or anybody on any forum we have tried. It does look like this is a dead end.

Now I ask you, can you shorten that to no more than a double integral with very few unrelated parameters? These are the types M can do, these are the only types anyone can do. Unless your supervisor is willing to give out a lot more information about this problem then maybe it is time to abandon this and look for something else. What have you got to lose? Please think about this.

zetafunc · 2017-02-21 21:53:17

Yes, I am trying to find out how to get a more reliable result. There are some ways to simplify the problem down to sums instead of integrals, or to use trigonometric functions instead of Bessel functions, but I don't know if that would help. I only have about 2 weeks to go until my project deadline, so I am willing to try anything until that point. Alternatively, we could instead get rid of the outer integral and instead figure out what the dependence on b,c is.

Last edited by zetafunc (2017-02-21 21:59:34)

bobbym · 2017-02-22 00:16:11

Alternatively, we could instead get rid of the outer integral and instead figure out what the dependence on b,c is.

How are you going to do that?

zetafunc · 2017-02-22 03:26:03

Well, the reason for the excessive number of variables is because we are expressing all of the vectors in terms of their components: for instance, we are taking k = (b,c). However, without loss of generality, we can set k = (r,0), r > 0. Then it suffices to investigate the integral

and plug in values of b to investigate how the value of the integral changes. However, I'm not sure if this simplification tells us something useful about the original problem -- though this is very easy to integrate with M and I have some results for how the integral changes as b varies between powers of 10.

bobbym · 2017-02-22 07:02:38

That is much better. It is a somewhat of a longshot but M has lots of stuff to find a relationship for a list of numbers.

zetafunc · 2017-02-22 07:34:56

I was planning on letting b = 10, b = 100, b = 1000, etc. and then looking at what the resulting integral evaluates to and dividing the ratios to determine the dependence the integral has on b. Is there a better way?

bobbym · 2017-02-22 14:02:46

Why powers of 10? Sometimes powers of 2 work. Is the data you are generating numerical?

zetafunc · 2017-02-22 21:37:05

Yes, it is purely numerical. My intuition says that the most significant differences will come from b large.

bobbym · 2017-02-22 22:00:40

You used the integral in post 12?

zetafunc · 2017-02-22 22:26:00

Yes, that one.

bobbym · 2017-02-22 22:36:01

What command in M?

zetafunc · 2017-02-24 09:13:14

The command is

NIntegrate[(BesselJ[1, x]*BesselJ[1, b - x])/(Abs[x]*
    Abs[b - x]), {x, -Infinity, Infinity}]

Some results:

b = 1: -0.307306
b = 2: 0.151302
b = 4: 0.29021
b = 8: -0.0672058
b = 16: 0.0372027
b = 32: -0.0154359
b = 64: -0.00547838
b = 128: -0.000120145
b = 256: 0.00056718
b = 512: 0.000169772
b = 1024: -0.0000548278
b = 2048: -0.0000149542

bobbym · 2017-02-24 12:09:52

Although it is oscillating around 0 it does look like as b -> infinity, that integral = 0.

Math Is Fun Forum

#1 2017-02-21 10:06:27

Another Bessel function integral

#2 2017-02-21 11:02:58

Re: Another Bessel function integral

#3 2017-02-21 11:06:01

Re: Another Bessel function integral

#4 2017-02-21 11:14:52

Re: Another Bessel function integral

#5 2017-02-21 11:22:45

Re: Another Bessel function integral

#6 2017-02-21 12:34:58

Re: Another Bessel function integral

#7 2017-02-21 12:44:36

Re: Another Bessel function integral

#8 2017-02-21 13:05:16

Re: Another Bessel function integral

#9 2017-02-21 15:43:15

Re: Another Bessel function integral

#10 2017-02-21 21:53:17

Re: Another Bessel function integral

#11 2017-02-22 00:16:11

Re: Another Bessel function integral

#12 2017-02-22 03:26:03

Re: Another Bessel function integral

#13 2017-02-22 07:02:38

Re: Another Bessel function integral

#14 2017-02-22 07:34:56

Re: Another Bessel function integral

#15 2017-02-22 14:02:46

Re: Another Bessel function integral

#16 2017-02-22 21:37:05

Re: Another Bessel function integral

#17 2017-02-22 22:00:40

Re: Another Bessel function integral

#18 2017-02-22 22:26:00

Re: Another Bessel function integral

#19 2017-02-22 22:36:01

Re: Another Bessel function integral

#20 2017-02-24 09:13:14

Re: Another Bessel function integral

#21 2017-02-24 12:09:52

Re: Another Bessel function integral

Board footer