Symbolic integration of Marcumq.
5 views (last 30 days)
Show older comments
Priyadarshi Mukherjee
on 13 Oct 2021
Commented: Priyadarshi Mukherjee
on 16 Oct 2021
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!
5 Comments
Walter Roberson
on 13 Oct 2021
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 ?
Accepted Answer
Walter Roberson
on 13 Oct 2021
Edited: Walter Roberson
on 13 Oct 2021
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);
8 Comments
Walter Roberson
on 15 Oct 2021
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)
inner_prod = prod(inner(2:end,:),1);
inner_prod(1:2)
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
More Answers (0)
See Also
Categories
Find more on Assumptions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!