Double integral over an anonymous function

1 visualizzazione (ultimi 30 giorni)
L'O.G.
L'O.G. il 6 Gen 2023
Risposto: Walter Roberson il 6 Gen 2023
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?

Risposte (1)

Walter Roberson
Walter Roberson il 6 Gen 2023
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.

Prodotti


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by