Error in double integration

I need to evaluate a double integration but showing some error

lam=532*10^-9;
z=100;
k=2*pi/lam;
omega=30;s=5;
r=linspace(0,100,100);
w0=0.002;m=1;rho=1;p=1;
E = 1./(w0^2) + (1i*k)./(2*z);
Con1=(1i./(2*lam*z))*exp((omega^2./2 + r.^2 + omega)./E).*exp((-1i.*k.*r.^2)./(2*z)).*(pi./E).*(1./(2*1i*sqrt(E)))^m;
Con2=(1i./(2*lam*z))*exp((omega^2./2 + r.^2 + omega)./conj(E)).*exp((-1i.*k.*r.^2)./(2*z)).*(pi./conj(E)).*(1./(2*1i*sqrt(conj(E)))).^m;
syms ph phi

for l=0:m
    E0 = Con1.*((exp(-r.*(cos(ph)+sin(ph))).*hermiteH(l,1i*(omega./2 - r*cos(ph))./sqrt(E)).*hermiteH(m-l,1i*(omega./2 - r*sin(ph))./sqrt(E))) - (exp(r.*(cos(ph)+sin(ph))).*hermiteH(l,-1i*(omega./2 + r*cos(ph))./sqrt(E)).*hermiteH(m-l,-1i*(omega./2 + r*sin(ph))./sqrt(E))));
    E0s = Con2.*((exp(-r.*(cos(phi)+sin(phi))).*hermiteH(l,1i*(omega./2 - r*cos(phi))./sqrt(conj(E))).*hermiteH(m-l,1i*(omega./2 - r*sin(phi))./sqrt(conj(E)))) - (exp(r.*(cos(phi)+sin(phi))).*hermiteH(l,-1i*(omega./2 + r.*cos(phi))./sqrt(conj(E))).*hermiteH(m-l,-1i*(omega./2 + r*sin(phi))./sqrt(conj(E)))));
    I = ((1i*p).^(m-l)).*nchoosek(m,l).*E0.*E0s.*exp(-1i.*s.*(ph-phi)).*exp(((-2.*r.^2)-(2.*r.^2.*cos(phi-ph)))./(rho.^2));
end
IQ = I;

ph_min = 0;
ph_max = 2*pi;
phi_min = 0;
phi_max = 2*pi;

fun = matlabFunction(IQ,'Vars',[ph,phi]);
result = integral2(fun, ph_min, ph_max, phi_min, phi_max);

Error using integral2Calc>integral2t/tensor
Integrand output size does not match the input size.

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Answer 1

[Torsten]

syms r ph phi
lam=532*10^-9;
z=100;
k=2*pi/lam;
omega=30;
w0=0.002;m=1;rho=1;p=1;
E = 1./(w0^2) + (1i*k)./(2*z);
Con1=(1i./(2*lam*z))*exp((omega^2./2 + r.^2 + omega)./E).*exp((-1i.*k.*r.^2)./(2*z)).*(pi./E).*(1./(2*1i*sqrt(E)))^m;
Con2=(1i./(2*lam*z))*exp((omega^2./2 + r.^2 + omega)./conj(E)).*exp((-1i.*k.*r.^2)./(2*z)).*(pi./conj(E)).*(1./(2*1i*sqrt(conj(E)))).^m;


for l=0:m
    E0 = Con1.*((exp(-r.*(cos(ph)+sin(ph))).*hermiteH(l,1i*(omega./2 - r*cos(ph))./sqrt(E)).*hermiteH(m-l,1i*(omega./2 - r*sin(ph))./sqrt(E))) - (exp(r.*(cos(ph)+sin(ph))).*hermiteH(l,-1i*(omega./2 + r*cos(ph))./sqrt(E)).*hermiteH(m-l,-1i*(omega./2 + r*sin(ph))./sqrt(E))));
    E0s = Con2.*((exp(-r.*(cos(phi)+sin(phi))).*hermiteH(l,1i*(omega./2 - r*cos(phi))./sqrt(conj(E))).*hermiteH(m-l,1i*(omega./2 - r*sin(phi))./sqrt(conj(E)))) - (exp(r.*(cos(phi)+sin(phi))).*hermiteH(l,-1i*(omega./2 + r.*cos(phi))./sqrt(conj(E))).*hermiteH(m-l,-1i*(omega./2 + r*sin(phi))./sqrt(conj(E)))));
    I = ((1i*p).^(m-l)).*nchoosek(m,l).*E0.*E0s.*exp(-1i.*l.*(ph-phi)).*exp(((-2.*r.^2)-(2.*r.^2.*cos(phi-ph)))./(rho.^2));
end
IQ = I;

ph_min = 0;
ph_max = 2*pi;
phi_min = 0;
phi_max = 2*pi;

fun = matlabFunction(IQ,'Vars',[ph,phi,r]);
r_array = linspace(0,10,100);
result = arrayfun(@(r)integral2(@(ph,phi)fun(ph,phi,r),ph_min, ph_max, phi_min, phi_max),r_array)

result = 
  Columns 1 through 10

  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 + 0.0000i

  Columns 11 through 20

  -0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 21 through 30

  -0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i

  Columns 31 through 40

   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i

  Columns 41 through 50

  -0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i

  Columns 51 through 60

   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0001 + 0.0001i   0.0000 - 0.0001i  -0.0001 - 0.0001i   0.0002 - 0.0001i  -0.0002 - 0.0002i  -0.0000 - 0.0003i  -0.0001 + 0.0004i  -0.0001 + 0.0005i

  Columns 61 through 70

  -0.0003 - 0.0007i  -0.0009 - 0.0001i   0.0003 - 0.0012i  -0.0016 + 0.0002i  -0.0017 - 0.0012i   0.0017 + 0.0021i   0.0027 + 0.0023i  -0.0047 - 0.0003i  -0.0032 + 0.0053i  -0.0068 - 0.0043i

  Columns 71 through 80

  -0.0013 + 0.0105i  -0.0099 + 0.0098i   0.0154 - 0.0098i   0.0166 - 0.0174i  -0.0018 + 0.0315i   0.0370 + 0.0187i  -0.0221 + 0.0498i   0.0713 - 0.0064i   0.0819 + 0.0464i  -0.0977 - 0.0761i

  Columns 81 through 90

  -0.1530 - 0.0558i   0.1956 - 0.0875i  -0.0220 - 0.2810i   0.3640 - 0.0717i  -0.3396 - 0.3510i  -0.1472 - 0.6259i   0.1278 + 0.8370i   0.5641 + 0.9619i  -1.4387 - 0.2965i  -0.8911 + 1.7179i

  Columns 91 through 100

  -2.4397 - 0.7423i   1.1402 + 3.1615i  -0.5514 + 4.3955i   0.7996 - 5.7850i  -2.3376 - 7.3367i   9.5097 + 3.5597i   7.1866 -11.3004i  16.7884 + 5.4952i  -8.9062 -21.5352i   0.3631 -30.7457i

Answer 2

[Surya]

Hi,

The error message suggests that the size of the output of the integrand function does not match the input size. This usually happens when there is a mismatch between the dimensions of the output and the dimensions expected by the integration function.

In your code, the integrand function fun has two input arguments ph and phi, and its output IQ is a complex-valued array. However, the dimensions of IQ are not consistent with the expected dimensions by integral2 function.

To fix the error, you may need to reshape the output of the integrand function IQ to match the expected dimensions. For example, if you expect a scalar output, you can modify your code as follows:

fun = matlabFunction(IQ,'Vars',[ph,phi]);
result = integral2(@(x,y) reshape(fun(x,y),[],1), ph_min, ph_max, phi_min, phi_max);

Here, the integral2 function expects a scalar output, so we use the reshape function to convert the output of the fun function to a column vector. The [] argument in reshape infers the size of the second dimension based on the size of the input, while the 1 specifies the size of the first dimension as 1.

Hope it helps.

[Athira T Das]

@Surya Modified my code, but still showing the same error