Symbolic integration of Marcumq.
6 ビュー (過去 30 日間)
古いコメントを表示
Priyadarshi Mukherjee
2021 年 10 月 13 日
コメント済み: Priyadarshi Mukherjee
2021 年 10 月 16 日
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 件のコメント
Walter Roberson
2021 年 10 月 13 日
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 ?
採用された回答
Walter Roberson
2021 年 10 月 13 日
編集済み: Walter Roberson
2021 年 10 月 13 日
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 件のコメント
Walter Roberson
2021 年 10 月 15 日
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
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Assumptions についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!