Symbolic integration of Marcumq.

Hello all.

I need to evaluate the integral of the product of a Bessel function and a Marcum Q function, which MATLAB cannot do numerically. So, the only option is to do it numerically in MATHEMATICA and use the value in my MATLAB code. But this appears to be very complex given the structure of my equation. Is there any way out?

Any suggestion will be helpful. Thanks!

[Walter Roberson]

MATLAB might possibly be able to work numerically ?

syms a b x Z
I0(Z) = besseli(0,Z)

I0(Z) = 

MQ(a,b) = int(x * exp((a^2-x^2)/2) * I0(a*x),x,b,inf)

MQ(a, b) = 

I2(Z) = besselj(2,Z)

I2(Z) = 

test_function = I2(Z) * MQ(a,b)

test_function = 

result = int(test_function, Z, 0, 1)

result = 

tr = subs(result,[a,b],[2,5])

tr = 

double(tr)

ans = 0.0048

[Priyadarshi Mukherjee]

Thanks a lot! I think this may work in my case now. I will try this approach now.

[Priyadarshi Mukherjee]

@Walter Roberson, I have the following problem to be specific.

We can see that the argument of Marcum Q is the parameter, which is to be integrated. This is creating all the difficulties.

[Walter Roberson]

I see a product over k = 2 to N, but I do not see k being used in the expression? In the form currently expressed it would be the same a taking [1-Q1(x,p)]^(N-1) instead of bothering with the product ?

[Priyadarshi Mukherjee]

Oh... I am sorry. In the process of simplifying my expression to state here, I omitted some index by mistake. Here is the thing:

Stuck in this for quite some time now.

On a different note, I was trying if the following helps:

Answer 1

[Walter Roberson]

Using the online tool here, execution can get through the calculation of "outer", but the integrals might take a while.

Could you fill out the Release (Version) field on the right? There are some approaches that became available in newer releases.

What is the order of magnitude of the sizes of the vector?

N = 10;
syms alpha [1 N] real
syms beta [1 N] real
syms gamma [1 N] real
syms delta [1 N] real
syms epsilon [1 N] real
syms a b x x_1 x_th Z real

I0(Z) = besseli(0,Z)
MQ1(a,b) = int(x * exp((a^2-x^2)/2) * I0(a*x),x,b,inf)

inner = repmat(1 - MQ1(delta(:) .* x_1, epsilon(:)), 1, N);
inner(1:N+1:end) = 1;
inner_prod = prod(inner(2:end,:),1);
COL = @(X) X(:);
outer = x_1 .* exp(-COL(beta(2:end)) .* x_1.^2) .* I0(COL(gamma(2:end))) .* inner_prod;
outer_int = int(outer, x_1, 0, x_th);
wanted_value = sum(COL(alpha(2:end)).*outer_int);

[Priyadarshi Mukherjee]

Thank you for your efforts.

N has the range 1:1000 and it cannot cross 1000 at any point of time.

But I could not understand what you meant by the following: "Could you fill out the Release (Version) field on the right? There are some approaches that became available in newer releases."

[Walter Roberson]

What matlab version are you using?

[Priyadarshi Mukherjee]

2018a

[Walter Roberson]

The below, with N = 10, executes in less than 10 seconds on my system. If you change the vpaintegral() to int() calls, then on my system it takes about 8 minutes.

The vpaintegral() and int() forms do not do the integrations -- the system is not one that MATLAB knows how to integrate to closed form, at least now when all the coefficients are symbolic.

To do numeric integration, you would have to put in numeric values for alpha, beta, gamma, delta, epsilon, a, b, x_1, and x_th (but leave x and Z symbolic.)

Do you have sample values that we could test with ?

tic
N = 10;
syms alpha [1 N] real
syms beta [1 N] real
syms gamma [1 N] real
syms delta [1 N] real
syms epsilon [1 N] real
syms a b x_1 x_th real
syms x Z real

I0(Z) = besseli(0,Z);
MQ1(a,b) = vpaintegral(x * exp((a^2-x^2)/2) * I0(a*x),x,b,inf);

inner = repmat(1 - MQ1(delta(:) .* x_1, epsilon(:)), 1, N);
inner(1:N+1:end) = 1;
inner_prod = prod(inner(2:end,:),1);
COL = @(X) X(:);
outer = x_1 .* exp(-COL(beta(2:end)) .* x_1.^2) .* I0(COL(gamma(2:end))) .* inner_prod;
outer_int = vpaintegral(outer, x_1, 0, x_th);
wanted_value = sum(COL(alpha(2:end)).*outer_int);

toc

[Walter Roberson]

If you are using Desktop, look to the right hand side of this posting, scroll down from the top.

If you are using mobile to access, you would need to scroll to near the bottom, below the contents, to see the field.

Notice the "Release" field. I have filled in R2018a for you.

R2018a tells me that you have vpaintegral() available to you, but that you do not have the much newer "hold" option supported by int() available to you.

[Priyadarshi Mukherjee]

Thanks a lot. I have some other values, from which I have to derive the values of these arrays.

1. In my MATLAB, "syms alpha [1 N] real" does not work, but the format "A = sym('a',[1 4])" works. So, we have to modify a bit.
2. Also, I think may be we have a problem in the following:

inner_prod = prod(inner(2:end,:),1);

COL = @(X) X(:);

outer = x_1 .* exp(-COL(beta(2:end)) .* x_1.^2) .* I0(COL(gamma(2:end))) .* inner_prod;

inner_prod will always be a 1*N all ones array.

3. We are most probably overlooking the "k not equal to i" condition.

4. I think we have to leave x_1 as a symbolic too, as it is the variable of integration.

But yes, I got the clue of how to approach the problem... thanks a lot (again)!

[Walter Roberson]

tic
N = 10;
alpha = sym('alpha', [1 N], 'real');
beta = sym('alpha', [1 N], 'real');
gamma = sym('alpha', [1 N], 'real');
delta = sym('alpha', [1 N], 'real');
epsilon = sym('alpha', [1 N], 'real');

syms a b x_th real
syms x x_1 Z real

I0(Z) = besseli(0,Z);
MQ1(a,b) = vpaintegral(x * exp((a^2-x^2)/2) * I0(a*x),x,b,inf);

inner = repmat(1 - MQ1(delta(:) .* x_1, epsilon(:)), 1, N);
inner(1:N+1:end) = 1;
inner(1:3,1:3)

ans = 

inner_prod = prod(inner(2:end,:),1);

inner_prod(1:2)

ans = 

Doesn't look like all 1's to me.

"We are most probably overlooking the "k not equal to i" condition"

Nope. Notice that I set the diagonal to 1. Multiplying A B C E is equivalent to multiplying A B C 1 E. With the diagonal set to 1, each column has effectively "crossed-out" the k == i value appropriate for its column.

Also, the bit about starting from 2 is why I index 2:end: it was easier to compute all the values including for index 1, and set the diagonal to 1's, multiply out, and ignore the first result, then to figure out the correct way to set the diagonal for a matrix that skipped the first row.

COL = @(X) X(:);
outer = x_1 .* exp(-COL(beta(2:end)) .* x_1.^2) .* I0(COL(gamma(2:end))) .* inner_prod;
outer_int = vpaintegral(outer, x_1, 0, x_th);
wanted_value = sum(COL(alpha(2:end)).*outer_int);

toc

Elapsed time is 2.884833 seconds.

[Priyadarshi Mukherjee]

Got it! I really appreciate your help. Thank you.