# Nested Numerical Integration with multivariable functions

11 views (last 30 days)
Scott Kaiser on 23 Aug 2023
Edited: Torsten on 23 Aug 2023
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 on 23 Aug 2023
Use integral3.

Torsten on 23 Aug 2023
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
##### 3 CommentsShow 1 older commentHide 1 older comment
Bruno Luong on 23 Aug 2023
Edited: Bruno Luong on 23 Aug 2023
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 on 23 Aug 2023
Edited: Torsten on 23 Aug 2023
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);

### Categories

Find more on Numerical Integration and Differentiation in Help Center and File Exchange

R2022a

### Community Treasure Hunt

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

Start Hunting!