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

Answer 1

[Alan Stevens]

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

[fatemeh kamali]

Thanks for your answer Alan Stevens! Yes i am sure about the definition of the functions.

[Alan Stevens]

Then it looks like the only sensible numerical upper limit for your integrals is k0.