How to calculate and plot ndefinite triple integral?

5 次查看（过去 30 天）
Hexe
Hexe 2023-4-10
评论： Hexe 2023-4-21
I have a triple indefinite integral (image attached).
Here respectively sx = sy = s*sin(a)/sqrt(2) and sz= s*cos(a). Parameter s=0.1 and parameter a changes from 0 to pi/2 – 10 points can be chosen [0 10 20 30 40 50 60 70 80 90]. Is it possible to solve such integral and to obtain the curve – plot(a,F)?
s=0.1;
a = 0:10:90;
fun = @(x,y,z) ((x.*z)./((x.^2+y.^2+z.^2))).*((2*pi)^(3/2))*exp(-(0.5.*sqrt(x.^2+y.^2+z.^2))).*exp(1i.*x*(s*sin(p)/sqrt(2))-2*((x.^2+y.^2+z.^2)+((z.^2)./((x.^2+y.^2+z.^2))))).*exp(1i.*y*(s*sin(p)/sqrt(2))-2*((x.^2+y.^2+z.^2)+((z.^2)./((x.^2+y.^2+z.^2))))).*exp(1i.*z*(s*cos(p))-2*((x.^2+y.^2+z.^2)+((z.^2)./((x.^2+y.^2+z.^2)))));
f3 = arrayfun(@(p)integral3(@(x,y,z)fun(x,y,z,p)),a);
plot(a,f3);
  3 个评论
显示 1更早的评论
Hexe
Hexe 2023-4-12
编辑：Hexe 2023-4-12
You are right. I forgot about the coefficient (2*pi)^(3/2) before exponent, but it does not matter much. The inportant thing is that in the second exponent there are 2 vectors: q and s. For the qx and qy sx=sy=s*sin(a)/sqrt(2) and for the qz sz=s*cos(a). Thus the code looks different than the written formula. Or the code for this case muct be written otherwise?
Thank you, I forgot about integration limits: [0, inf, 0, 2*pi, 0, pi].
Torsten
Torsten 2023-4-12
I forgot about the coefficient (2*pi)^(3/2) before exponent, but it does not matter much.
There are many more differences.
In your formula:
exp(-0.5.*(x.^2+y.^2+z.^2))
In your code:
exp(-(0.5.*sqrt(x.^2+y.^2+z.^2)))
In your formula:
exp(1i.*x*(s*sin(p)/sqrt(2))+1i.*y*(s*sin(p)/sqrt(2))+1i*z.*(s*cos(p))-2*(x.^2+y.^2+z.^2+z.^2./(x.^2+y.^2+z.^2)))
In your code:
exp(1i.*x*(s*sin(p)/sqrt(2))-2*((x.^2+y.^2+z.^2)+((z.^2)./((x.^2+y.^2+z.^2))))).*exp(1i.*y*(s*sin(p)/sqrt(2))-2*((x.^2+y.^2+z.^2)+((z.^2)./((x.^2+y.^2+z.^2))))).*exp(1i.*z*(s*cos(p))-2*((x.^2+y.^2+z.^2)+((z.^2)./((x.^2+y.^2+z.^2)))))

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-4-12
Ran in:
s = 0.1;
a = 0:5:360;
a = a*pi/180;
fun = @(x,y,z,p) x.*z./(x.^2+y.^2+z.^2).*exp(-0.5*(x.^2+y.^2+z.^2)).*exp(1i*x*(s*sin(p)/sqrt(2))+1i*y*(s*sin(p)/sqrt(2))+1i*z*(s*cos(p))-2*(x.^2+y.^2+z.^2+z.^2./(x.^2+y.^2+z.^2)));
f3 = (2*pi)^1.5*arrayfun(@(p)integral3(@(x,y,z)fun(x,y,z,p),0,Inf,0,2*pi,0,pi),a);
figure(1)
plot(a,real(f3))
figure(2)
plot(a,imag(f3))
  8 个评论
显示 6更早的评论
Torsten
Torsten 2023-4-19
编辑：Torsten 2023-4-19
Ran in:
Why do you replace s by k and not by m in your code ?
And if you loop over the elements of a, why do you use the arrayfun ? Arrayfun computes the values for f3 for the complete vector a over and over again. I can understand that your code takes a while to finish.
Since the results for f3 are complex-valued, you can only apply surf on abs(f3) or imag(f3) or real(f3), but not f3 itself.
n = 1;
t = 1;
r = 1;
S = 1:0.5:5;
P = 0:10:180;
P = P*pi/180;
for i = 1:numel(S)
s = S(i);
for j = 1:numel(P)
p = P(j);
fun = @(x,y,z) x.*z./(x.^2+y.^2+z.^2).*exp(-0.5*(x.^2+y.^2+z.^2)).*exp(1i*x*(s*sin(p)/sqrt(2))+1i*y*(s*sin(p)/sqrt(2))+1i*z* (s*cos(p))-2*(x.^2+y.^2+z.^2+z.^2./(x.^2+y.^2+z.^2)));
f3(i,j) = (2*pi)^1.5*integral3(fun,0,Inf,0,2*pi,0,pi);
end
end
f3
f3 =
Columns 1 through 10 0.2759 + 0.0978i 0.2654 + 0.1215i 0.2546 + 0.1408i 0.2450 + 0.1555i 0.2375 + 0.1655i 0.2330 + 0.1710i 0.2320 + 0.1722i 0.2346 + 0.1689i 0.2405 + 0.1613i 0.2490 + 0.1491i 0.2481 + 0.1379i 0.2262 + 0.1685i 0.2041 + 0.1917i 0.1847 + 0.2080i 0.1700 + 0.2183i 0.1613 + 0.2237i 0.1594 + 0.2247i 0.1644 + 0.2215i 0.1759 + 0.2138i 0.1929 + 0.2008i 0.2134 + 0.1689i 0.1782 + 0.2014i 0.1439 + 0.2231i 0.1145 + 0.2360i 0.0928 + 0.2427i 0.0802 + 0.2455i 0.0776 + 0.2459i 0.0849 + 0.2442i 0.1017 + 0.2396i 0.1268 + 0.2302i 0.1750 + 0.1897i 0.1269 + 0.2193i 0.0816 + 0.2345i 0.0444 + 0.2395i 0.0179 + 0.2392i 0.0031 + 0.2375i 0.0001 + 0.2369i 0.0087 + 0.2377i 0.0288 + 0.2387i 0.0599 + 0.2368i 0.1360 + 0.2006i 0.0769 + 0.2231i 0.0239 + 0.2278i -0.0173 + 0.2221i -0.0451 + 0.2130i -0.0598 + 0.2061i -0.0627 + 0.2044i -0.0540 + 0.2083i -0.0334 + 0.2160i -0.0002 + 0.2235i 0.0991 + 0.2028i 0.0321 + 0.2151i -0.0244 + 0.2074i -0.0651 + 0.1900i -0.0903 + 0.1722i -0.1028 + 0.1605i -0.1050 + 0.1579i -0.0976 + 0.1646i -0.0794 + 0.1787i -0.0481 + 0.1958i 0.0660 + 0.1979i -0.0053 + 0.1987i -0.0608 + 0.1782i -0.0967 + 0.1502i -0.1163 + 0.1255i -0.1247 + 0.1104i -0.1259 + 0.1072i -0.1208 + 0.1159i -0.1073 + 0.1349i -0.0814 + 0.1600i 0.0379 + 0.1880i -0.0342 + 0.1771i -0.0849 + 0.1451i -0.1130 + 0.1090i -0.1249 + 0.0802i -0.1284 + 0.0635i -0.1285 + 0.0602i -0.1263 + 0.0699i -0.1188 + 0.0914i -0.1004 + 0.1218i 0.0149 + 0.1748i -0.0549 + 0.1531i -0.0981 + 0.1121i -0.1165 + 0.0712i -0.1201 + 0.0412i -0.1188 + 0.0250i -0.1180 + 0.0220i -0.1186 + 0.0315i -0.1172 + 0.0532i -0.1072 + 0.0857i Columns 11 through 19 0.2593 + 0.1324i 0.2701 + 0.1111i 0.2801 + 0.0856i 0.2880 + 0.0566i 0.2929 + 0.0252i 0.2940 - 0.0074i 0.2912 - 0.0397i 0.2849 - 0.0702i 0.2759 - 0.0978i 0.2136 + 0.1816i 0.2357 + 0.1553i 0.2566 + 0.1218i 0.2736 + 0.0818i 0.2841 + 0.0369i 0.2866 - 0.0104i 0.2807 - 0.0571i 0.2672 - 0.1003i 0.2481 - 0.1379i 0.1583 + 0.2137i 0.1930 + 0.1877i 0.2267 + 0.1509i 0.2546 + 0.1034i 0.2724 + 0.0476i 0.2768 - 0.0125i 0.2670 - 0.0718i 0.2446 - 0.1251i 0.2134 - 0.1689i 0.1001 + 0.2277i 0.1462 + 0.2071i 0.1925 + 0.1718i 0.2322 + 0.1209i 0.2581 + 0.0570i 0.2649 - 0.0137i 0.2510 - 0.0833i 0.2189 - 0.1438i 0.1750 - 0.1897i 0.0450 + 0.2250i 0.0994 + 0.2139i 0.1566 + 0.1843i 0.2075 + 0.1338i 0.2418 + 0.0651i 0.2514 - 0.0140i 0.2335 - 0.0916i 0.1920 - 0.1564i 0.1360 - 0.2006i -0.0024 + 0.2090i 0.0562 + 0.2094i 0.1212 + 0.1888i 0.1817 + 0.1421i 0.2241 + 0.0715i 0.2368 - 0.0133i 0.2154 - 0.0967i 0.1652 - 0.1631i 0.0991 - 0.2028i -0.0395 + 0.1839i 0.0189 + 0.1963i 0.0882 + 0.1862i 0.1560 + 0.1461i 0.2055 + 0.0763i 0.2215 - 0.0118i 0.1975 - 0.0989i 0.1400 - 0.1649i 0.0660 - 0.1979i -0.0654 + 0.1541i -0.0110 + 0.1773i 0.0589 + 0.1781i 0.1313 + 0.1461i 0.1867 + 0.0796i 0.2060 - 0.0096i 0.1803 - 0.0986i 0.1172 - 0.1627i 0.0379 - 0.1880i -0.0810 + 0.1233i -0.0332 + 0.1550i 0.0340 + 0.1660i 0.1082 + 0.1428i 0.1680 + 0.0814i 0.1906 - 0.0070i 0.1642 - 0.0964i 0.0971 - 0.1576i 0.0149 - 0.1748i
Hexe
Hexe 2023-4-21
Thank you very much. Your notes helped me to build the necessary surface.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Graphics Performance 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by