Nested Numerical Integration with multivariable functions

8 次查看(过去 30 天)
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
Torsten 2023-8-23
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 个评论
Bruno Luong
Bruno Luong 2023-8-23
编辑:Bruno Luong 2023-8-23
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
Torsten 2023-8-23
编辑:Torsten 2023-8-23
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);

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by