Matlab Numerical integral issue
显示 更早的评论
Hi everyone.
I'm trying to solve the integral in the code below:
clc
clear all
%Enviroment Properties
f=15e9;
c=2.99792458e8;
lambda=c/f;
mu0=4*pi*1e-7;
eps0=1.0/c^2/mu0;
eta0=sqrt(mu0/eps0);
epsr=5;
omega=2*pi*f;
k0=omega*sqrt(mu0*eps0);
k2=omega*sqrt(mu0*eps0*epsr);
%square patch array properties
l=1.8e-3;
alpha_Ett=1.02*(l^3);
p=2e-3;
N=1/(p^2);
r=0.6956*p;
chieexx=(N*alpha_Ett)/(1-(N*alpha_Ett/(4*r)));
%source and observation points
z0=3e-3;
z=2e-3;
y=0;
for x=5e-6:1e-5:5e-3
pho=sqrt(x.^2+y.^2);
f1 = (@(kp) (kp.^3./(sqrt(k0^2-kp.^2))).* besselj(0,kp.*pho).* (((2*1i*(sqrt(k0^2-kp.^2)).*(k2^2))+ ((sqrt(k2^2-kp.^2)).*(((sqrt(k0^2-kp.^2))*(k2^2).*chieexx)+ ((k0^2)*(-2i+(sqrt(k0^2-kp.^2)).*chieexx)))))./((2*1i*(sqrt(k0^2-kp.^2)).*(k2^2))+ ((sqrt(k2^2-kp.^2)).*(((sqrt(k0^2-kp.^2)).*(k2^2).*chieexx)+ ((k0^2).*(2i+(sqrt(k0^2-kp.^2)).*chieexx)))))) .* exp(1i.*(sqrt(k0^2-kp.^2)).*(z+z0)) );
f2= (@(kp) (kp.^3./(sqrt(k0^2-kp.^2))).* besselj(0,kp.*pho).* exp(1i.*(sqrt(k0^2-kp.^2)).*(-z+z0)) );
I1= quadgk(f1,0,inf,'RelTol',1e-8);
I2= quadgk(f2,0,inf,'RelTol',1e-8);
E= (eta0/(4*pi*k0)).*(abs(I1+I2));
semilogy(x,abs(E),'Or')
hold on
end
set(gca,'FontSize',14,'LineWidth',3)
xlabel('x [m]');ylabel('|Ez| [V/m]')
But i get the following warning :
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.2e+03.
The integral may not exist, or it may be difficult to approximate numerically. Increase MaxIntervalCount to
1132 to enable QUADGK to continue for another iteration.
Could you please help me to improve my code.
Thanks in advance
回答（1 个）
Alan Stevens
2021-3-13
编辑：Alan Stevens
2021-3-13
Your problem occurs when the upper limit of the integration exceeds k0. Are you sure the functions are defined correctly for that situation?
If the upper limit is taken as k0, the following code
%Enviroment Properties
f=15e9;
c=2.99792458e8;
lambda=c/f;
mu0=4*pi*1e-7;
eps0=1.0/c^2/mu0;
eta0=sqrt(mu0/eps0);
epsr=5;
omega=2*pi*f;
k0=omega*sqrt(mu0*eps0);
k2=omega*sqrt(mu0*eps0*epsr);
%square patch array properties
l=1.8e-3;
alpha_Ett=1.02*(l^3);
p=2e-3;
N=1/(p^2);
r=0.6956*p;
chieexx=(N*alpha_Ett)/(1-(N*alpha_Ett/(4*r)));
%source and observation points
z0=3e-3;
z=2e-3;
y=0;
%x = 5e-6:1e-4:5e-3;
x = logspace(-5.301,-2.301);
lo = 0; hi = k0; % Integration limits
for i = 1:numel(x)
pho=sqrt(x(i).^2+y.^2);
f1 = (@(kp) (kp.^3./(sqrt(k0^2-kp.^2))).* besselj(0,kp.*pho).* (((2*1i*(sqrt(k0^2-kp.^2)).*(k2^2))+ ((sqrt(k2^2-kp.^2)).*(((sqrt(k0^2-kp.^2))*(k2^2).*chieexx)+ ((k0^2)*(-2i+(sqrt(k0^2-kp.^2)).*chieexx)))))./((2*1i*(sqrt(k0^2-kp.^2)).*(k2^2))+ ((sqrt(k2^2-kp.^2)).*(((sqrt(k0^2-kp.^2)).*(k2^2).*chieexx)+ ((k0^2).*(2i+(sqrt(k0^2-kp.^2)).*chieexx)))))) .* exp(1i.*(sqrt(k0^2-kp.^2)).*(z+z0)) );
f2= (@(kp) (kp.^3./(sqrt(k0^2-kp.^2))).* besselj(0,kp.*pho).* exp(1i.*(sqrt(k0^2-kp.^2)).*(-z+z0)) );
I1= quadgk(f1,lo,hi);
I2= quadgk(f2,lo,hi);
E(i)= (eta0/(4*pi*k0)).*(abs(I1+I2));
hold on
end
semilogy(x,abs(E),'Or')
set(gca,'FontSize',14,'LineWidth',3)
xlabel('x [m]');ylabel('|Ez| [V/m]')
produces this result
2 个评论
fatemeh kamali
2021-3-15
Alan Stevens
2021-3-15
Then it looks like the only sensible numerical upper limit for your integrals is k0.
类别
在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息
产品
2021-3-12
2021-3-15
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
