필터 지우기
필터 지우기

Double integral over an anonymous function

조회 수: 4 (최근 30일)
L'O.G.
L'O.G. 2023년 1월 6일
답변: Walter Roberson 2023년 1월 6일
I am trying to do the following double integral:
where
I attempt this by
q = 0.01:0.01:50; % range is somewhat arbitrary
fun = @(q) integral2(@(u,sigma) (((sin(q.*(B/2)*sqrt(1-sigma.^2).*cos(pi*u./2))./(q.*(B/2)*sqrt(1-sigma.^2).*cos(pi*u./2))).*(sin(q.*(B/2)*sqrt(1-sigma.^2).*a*sin(pi*u./2))./(q.*(B/2)*sqrt(1-sigma.^2).*a*sin(pi*u./2)))).^2).*(sin(q.*B*c*sigma./2)./(q.*B*c*sigma./2)).^2,0,1,0,1);
pq = arrayfun(@(q)fun(q),q);
I imagine I made a mistake, especially regarding the element wise operations, but I don't know how to write this in a cleaner way because that breaks up the anonymous function. Can somebody help me clean this up?

답변 (1개)

Walter Roberson
Walter Roberson 2023년 1월 6일
We could probably produce a more compact function if we knew the sign() of the various constants
syms a B c mu q sigma x real
assume(a ~= 0)
assume(B ~= 0)
assume(c ~= 0)
assume(q > 0) %an important restriction that avoids division by 0
assume(sigma ~= 0)
Pi = sym(pi);
S(x) = sin(x)/x
S(x) = 
inner_phi = (S(mu/2*cos(Pi/2*mu)) * S(mu*a/2 * sin(Pi/2*mu)))
inner_phi = 
phi(mu,a) = int(inner_phi^2, mu, 0, 1)
phi(mu, a) = 
Pmu_inner = phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2
Pmu_inner = 
Pmu(mu) = int(Pmu_inner, sigma, 0, 1)
Pmu(mu) = 
P(q) = Pmu(q*B)
P(q) = 
Pfun = matlabFunction(P, 'vars', [q, a, B, c, sigma])
Pfun = function_handle with value:
@(q,a,B,c,sigma)integral(@(sigma)1.0./B.^2.*1.0./c.^2.*1.0./q.^2.*1.0./sigma.^2.*sin((B.*c.*q.*sigma)./2.0).^2.*integral(@(mu)1.0./a.^2.*1.0./mu.^4.*sin((a.*mu.*sin((mu.*pi)./2.0))./2.0).^2.*1.0./cos((mu.*pi)./2.0).^2.*1.0./sin((mu.*pi)./2.0).^2.*sin((mu.*cos((mu.*pi)./2.0))./2.0).^2.*1.6e+1,0.0,1.0).*4.0,0.0,1.0)
Is this going to be usable in practice?
Possibly not. You have S() of expressions involving sin(mu) where mu is q*B and q ranges 0 to 50. Unless B is in the range 0 to 1/50 or so, then you are going to be taking sin() of an expression that ranges further than pi/2 so you are going to (at least in theory) be getting a 0 from that part. But it is an argument to S which is sin(x)/x and the x part is the sin(q*B*pi/2) and if that inner sin() can go to 0, then you are going to have division by 0, which is going to give you problems during numeric evaluation.
The sinc() function itself does not have problems with division by 0 because mathematically it is not just plain sin(x)/x but rather limit(sin(x+delta))/(x+delta) as delta->0 and the limit() goes to 1 there because near 0 sin(x) can be closely approximated by x, and limit (x+delta)/(x+delta) goes to 1.

카테고리

Help CenterFile Exchange에서 Logical에 대해 자세히 알아보기

제품


릴리스

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by