Integrating function defined as an integral with integral2

Hi,

I have a problem with integrating a function with integral2. The funtion I want to integrate is defined as the product of two functions, one of which is already defined as a double integral.

Some code to get in context :

%% Init
clc

%% Parameters
index = 1.4581;
w_eta = 20e-6;
w_neta = 20e-6;
w_f = 7.5e-6;
lbda = 1.55e-6;
k = 2*pi/lbda;
z1 = 100e-6;
Rx_lens = 20e-6;
Ry_lens = 20e-6;
z2 = 100e-6;

%% Functions
fun_U_LD = @(eta,neta) exp(-((eta./w_eta).^2 + (neta./w_neta).^2));

fun_U_fiber = @(u,v) exp(-(u.^2 + v.^2)./w_f.^2);

fun_conj_U_fiber = @(u,v) conj(exp(-(u.^2 + v.^2)./w_f.^2));

fun_U_lens = @(x,y) exp(1i*k*z1)/(1i*lbda*z1) ...
    .* integral2(@(eta,neta) fun_U_LD(eta,neta).* exp(1i*k/(2*z1) ...
    * ((x - eta).^2 + (y - neta).^2)), -Inf, Inf, -Inf, Inf);

fun_Delta = @(x,y) x.^2./(Rx_lens + sqrt(Rx_lens.^2 - x.^2)) ...
    + y.^2./(Ry_lens + sqrt(Ry_lens.^2 - y.^2));

fun_t_lens = @(x,y) exp(-1i*k*(index-1) ...
    .* (x.^2./(Rx_lens + sqrt(Rx_lens.^2 - x.^2)) ...
    + y.^2./(Ry_lens + sqrt(Ry_lens.^2 - y.^2))));

fun_U_t = @(x,y) fun_U_lens(x,y).*fun_t_lens(x, y);

integral2(fun_U_t, -Inf, Inf, -Inf, Inf)

And then these are the errors (a lot) I get :

Matrix dimensions must agree.

Error in coupl_eff_chao_hu>@(eta,neta)fun_U_LD(eta,neta).*exp(1i*k/(2*z1)*((x-eta).^2+(y-neta).^2))

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
    @(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc/iterateScalarValued (line 323)
                fx = FUN(t).*w;

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in
integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x))
(line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

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

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

Error in
coupl_eff_chao_hu>@(x,y)exp(1i*k*z1)/(1i*lbda*z1).*integral2(@(eta,neta)fun_U_LD(eta,neta).*exp(1i*k/(2*z1)*((x-eta).^2+(y-neta).^2)),-Inf,Inf,-Inf,Inf)

Error in @(x,y)fun_U_lens(x,y).*fun_t_lens(x,y)

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
    @(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in
integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x))
(line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

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

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

The main error seems to be with matrix dimensions, which is impossible to specifiy.

In another similar code I had with integral (and not integral2) I could calculate integrals of functions defined as integrals using the 'ArrayValued' option, but it doesn't exist for integral2.

I realize the definitions are complicated, and maybe I am trying something impossible but I am merely writing equations from a scientific paper into Matlab, and I would like to know if there is something I am doing wrong here. Please any advice, or even saying that this is impossible, would be appreciated, as I could move on to try another method.

Thanks

[Star Strider]

That is somewhat difficult to follow.

My only suggestion is to experiment with doing the double integration as ‘nested’ integrations using integral with the 'ArrayValued' option selected. That has worked for me in the past, however it may not be generally applicable (so likely the reason that it is not a part of the integral2 options), so I cannot say if it is applicable to your project.

[David Goodmanson]

HI Xavier,

it appears that all these functions are separable, e.g.

fun_t_lens = @(x,y) exp(-1i*k*(index-1) ...
    .* (x.^2./(Rx_lens + sqrt(Rx_lens.^2 - x.^2)) ...
    + y.^2./(Ry_lens + sqrt(Ry_lens.^2 - y.^2))));

and consequently

fun_t_lens = fun_t_lens_x * fun_t_lens_y
% where
fun_t_lens_x = @(x) exp(-1i*k*(index-1) .* (x.^2./(Rx_lens + sqrt(Rx_lens.^2 - x.^2))
fun_t_lens_y = @(y) exp(-1i*k*(index-1) .* (y.^2./(Ry_lens + sqrt(Ry_lens.^2 - y.^2))

moreover, eta is always associated with x and neta is always associated with y. If you separate the functions of eta and neta in the same way, it looks like at each stage you could reduce this to

integral(f(x),inf,inf) * integral(f(y),inf,inf)
instead of
integral2(f(x,y),inf,inf,inf,inf)

and similarly for

integral2(g(eta,neta),inf,inf,inf,inf)

There will still probably be issues with inserting integrals as arguments of other integrals, but at least it looks like the integral2 stuff can go away. Also, fun_U_lens can be integrated analytically by completing the square in the exponent.

[Xavier Insou]

Hello David and Star, thank you for your answers.

I will try to do the separate version with integral today if I get the time, and get back to you with the results.

Ultimately the goal is to calcultate the coupling efficiency between two optical beams, so with a known configuration I should know pretty quick if its working the way I want it to or not.

[Xavier Insou]

Hello, it didn't work with single integrals, I don't get errors, but infinite integrals send out warnings and take way too long. Plus the result isn't what expected.

Anyway thank you for your kind help. I am moving on to another method.