How to calculate and plot ndefinite triple integral?

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);

[Torsten]

The function in the image and in your code differ substantially. Further you didn't specify the limits of integration.

[Hexe]

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]

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)))))

Answer 1

[Torsten]

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))

[Hexe]

You are right, there are many differences that I did not see after writing code many times in different ways. I did not now how to write code for vectors in exponents. Thank you very much.

[Hexe]

I made a plot of the curve with coefficients and now it works right. But I also wanted to ask - is it possible to present the plot as a surface?

[Torsten]

But I also wanted to ask - is it possible to present the plot as a surface?

How ? f3 is a function of a single variable a. Where do you see a surface here ?

[Hexe]

Actually, s is a variable too, but due to the difficult function I wanted at first to obtain curves with fixed s - curve for s = 0.1, another one for for s =1 etc.

[Torsten]

Then loop over the values of s and a, save the results for f3 in a matrix f3(i,j) where i stands for a(i) and j stands for s(j) and make a surface plot:

surf(s,a,f3)

What's the problem ?

[Hexe]

I tried to make it, and Matlab calculates, but it too long and does not give a result. What is wrong with it?

n = 1;

t = 1;

r = 1;

s = 1:0.5:10;

a = 0:10:360;

a = a*pi/180;

for k = 1:length(s)

m = s(k);

for j = 1:length(a)

p = a(j);

fun = @(x,y,z,p) (x.*z./(x.^2+y.^2+z.^2)).*exp(-(0.5*(r^2).*(x.^2+y.^2+z.^2))).*exp((1i*k*sin(p)*(2^(-1/2))*x+1i*k*sin(p)*(2^(-1/2))*y+1i*k*cos(p)*z)-2*t*(x.^2+y.^2+z.^2+n*z.^2./(x.^2+y.^2+z.^2)));

f3(j,:) = arrayfun(@(p)integral3(@(x,y,z)fun(x,y,z,p),0,Inf,0,2*pi,0,pi),a);

Cor(k,j) = (-n*((1+4*t/(r^2))^(3/2))*2*sqrt(2*n*t)/(sqrt(pi)*erf(sqrt(2*n*t)))*(r^3)*(2*pi)^(3/2)/(2*pi)^3).*f3(k,j);

end

end

surf(a,s,Cor);

[Torsten]

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]

Thank you very much. Your notes helped me to build the necessary surface.