Nested Numerical Integration with multivariable functions

I am attempting to evaluate a nested integral in matlab. The equation comes from a paper I am reading and is a weighting function used within another integral. The weighting function I am trying to numerically calculate is the following:

The variables D,d, and L are constant. I tried symbolic integration, beginning with the inner integral, but the code does not perform the integration. Only returning the input with the integration bounds. I started looking into numerical integration but the trapz() function doesnt work well with symbols. My code is shown below. And of course this code doesnt work, but the other methods I have tried dont either. Is there a way to alter this code to work for me? Any suggestions?

clear;clc;

d=.01; %1/e patch diameter;

D=.35; %Aperture Diameter

L=7e3; %Propagation distance;

N=30; %number of terms in the numerical integration

delta_theta=.1; %[rad] angular separation between 2 patches. At L=7000m .1 rad is 70cm. (Like the width of a window)

syms x u theta z delta

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%Calculate the 1st term

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Define the coefficients

coeff=-11.64*(16/pi)^2*D^(-1/3);

%Define the 1st integrand

f1=x*exp(-2*x^2);

%Define the 3rd integrand

f3=((u*acos(u))-u^2*(3-2*u^2)*sqrt(1-u^2))*(u^2*(1-z/L)^2+(z*d*x/(D*L))^2+2*u*(1-z/L)*...

(z*d*x/(D*L))*cos(theta))^(5/6);

u_values=linspace(0,1,N);

inner=zeros(1, N);

for i=1:N

inner_integral=subs(f3,u,u_values(i));

inner(i)=inner_integral;

end

inner=trapz(u_values,inner);

%Evaluate the middle integral numerically

theta_values=linspace(0,2*pi,N);

middle=zeros(1, N);

for j=1:N

middle_integral=subs(inner(j),theta,theta_values(j));

middle(i)=middle_integral;

end

middle=trapz(theta_values,middle);

%Evaluate the outer integral numerically

x_values=linspace(0,10000,N);

outer=zeros(1, N);

for j=1:N

outer_integral=subs(f1*inner(j),x,x_values(j));

outer(i)=outer_integral;

end

outer=trapz(x_values,outer);

first_term=coeff*outer;

[Torsten]

Use integral3.

Answer 1

[Torsten]

d=.01; %1/e patch diameter;
D=.35;  %Aperture Diameter
L=7e3; %Propagation distance;
coeff=-11.64*(16/pi)^2*D^(-1/3);
fun = @(x,u,theta,z)x.*exp(-2*x.^2).*(u.*acos(u)-u.^2.*(3-2*u.^2).*sqrt(1-u.^2)).*(u.^2.*(1-z/L).^2+(d*z*x/(D*L)).^2+2*u.*(1-z/L).*...
                    (z*d*x/(D*L)).*cos(theta)).^(5/6);
z = 0:50;
fz = coeff*arrayfun(@(z)integral3(@(x,u,theta)fun(x,u,theta,z),0,Inf,0,1,0,2*pi),z)

fz = 1×51

17.2324 17.2283 17.2242 17.2201 17.2160 17.2119 17.2078 17.2037 17.1996 17.1955 17.1914 17.1873 17.1832 17.1791 17.1750 17.1709 17.1668 17.1627 17.1586 17.1545 17.1504 17.1464 17.1423 17.1382 17.1341 17.1300 17.1259 17.1218 17.1177 17.1136

plot(z,fz)
grid on

[Scott Kaiser]

Thank you for this! it seems efficient, and appears to work well. As a followup there is a second half to this weighting function that is similar to the first half given by "fun".

It is:

fun2 = @(delta,x,u,theta,z)x.*exp(-2*x.^2)...
    .*(u.*acos(u)-u.^2*(3-2*u.^2)*sqrt(1-u.^2))...
   .*(u.^2.*(1-z/L).^2+(z*d/(D*L))^2*(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))...
   +2*u.*(1-z/L)*(z*d/(D*L))*cos(theta)*sqrt(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))).^(5/6);

This is a function of 4 variables and is 4 nested integrals. I tried implementing your code with this, but I think integral3 works for only 3 variables. Is this true, or is there a work-around?

clear z;
coeff2=(5.82/pi)*(16/pi)^2*D^(-1/3);
fun2 = @(delta,x,u,theta,z)x.*exp(-2*x.^2)...
    .*(u.*acos(u)-u.^2*(3-2*u.^2)*sqrt(1-u.^2))...
   .*(u.^2.*(1-z/L).^2+(z*d/(D*L))^2*(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))...
   +2*u.*(1-z/L)*(z*d/(D*L))*cos(theta)*sqrt(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))).^(5/6);

z = linspace(0,L,N);
fz2 = coeff2*arrayfun(@(z)integral3(@(delta,x,u,theta)fun2(delta,x,u,theta,z),0,2*pi,0,Inf,0,1,0,2*pi),z);

Do you have any suggestions for working through this function "fun2"??

[Bruno Luong]

For 4D integration you can call

• integral over integral3 or
• integral3 over integral or
• integral2 over integral2
• trapz over integral3 (for example you can call trapz over to compute )
• etc...

[Torsten]

coeff2=(5.82/pi)*(16/pi)^2*D^(-1/3);
Z = linspace(0,L,N);
for i = 1:numel(Z)
    z = Z(i);
    fun2 = @(delta,x,u,theta)x.*exp(-2*x.^2)...
        .*(u.*acos(u)-u.^2*(3-2*u.^2)*sqrt(1-u.^2))...
        .*(u.^2.*(1-z/L).^2+(z*d/(D*L))^2*(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))...
        +2*u.*(1-z/L)*(z*d/(D*L))*cos(theta)*sqrt(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))).^(5/6);
    fz2(i) = coeff2*integral(@(x)integral3(@(delta,u,theta)fun2(delta,x,u,theta),0,2*pi,0,1,0,2*pi),0,Inf,'ArrayValued',true);
end

or maybe (if you think you can decipher in 2 weeks what you have done in the last line):

coeff2=(5.82/pi)*(16/pi)^2*D^(-1/3);
z = linspace(0,L,N);
fun2 = @(delta,x,u,theta,z)x.*exp(-2*x.^2)...
        .*(u.*acos(u)-u.^2*(3-2*u.^2)*sqrt(1-u.^2))...
        .*(u.^2.*(1-z/L).^2+(z*d/(D*L))^2*(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))...
        +2*u.*(1-z/L)*(z*d/(D*L))*cos(theta)*sqrt(x.^2+(L/d*delta_theta).^2-2*x*L/d*delta_theta.*cos(delta))).^(5/6);
fz2 = coeff2*arrayfun(@(z)integral(@(x)integral3(@(delta,u,theta)fun2(delta,x,u,theta,z),0,2*pi,0,1,0,2*pi),0,Inf,'ArrayValued',true),z);