triple nested integral with error "Arrays have incompatible sizes for this operation."

Hi! I have a code in which there are three integrals, nested inside one another. These are the functions:

function integral_z = integral_z(z0, a, y, P, A, omega0, to)
    fe = f_e(z0, P);
    insideF0 = 2.*a.*(1+7.*(a.^2)).*(2*pi*P.nu/omega0).*(to.^(3/2)).*(z0.^(-1) + A.*(8.*to).^(-1/2).*a.^(-1).*y.^(-2).*(1-y.^(16/9))).^2;
    integral_z = fe.*F0(insideF0);
end

function Integrand = Integrand(a, y, P, A, omega0, to)
    z0_min = ((P.p-2)./(P.p-1))./(1-P.geratio.^(2-P.p));
    z0_max = P.geratio.*z0_min;
    % calling first integral
    Integrand = ((a.^3).*((1+7.*(a.^2)).^(-2))).*integral(@(z0)integral_z(z0, a, y, P, A, omega0, to), z0_min, z0_max);
end

function IntegrandY = IntegrandY(y, P, A, omega0, to)
    % calling second integral
    IntegrandY = y.^3.*integral(@(a)Integrand(a, y, P, A, omega0, to), 0, 1);
end

and I call IntegrandY as follows:

% calling third integral
integral(@(y)IntegrandY(y, P, A, omega0, to),0,1);

But I get the following error:

Arrays have incompatible sizes for this operation.

Error in F_BM_Analytic_Cooling/integral_z (line 49)
    insideF0 = 2.*a.*(1+7.*(a.^2)).*(2*pi*P.nu/omega0).*(to.^(3/2)).*(z0.^(-1) + A.*(8.*to).^(-1/2).*a.^(-1).*y.^(-2).*(1-y.^(16/9))).^2;

Error in F_BM_Analytic_Cooling>@(z0)integral_z(z0,a,y,P,A,omega0,to) (line 56)
    Integrand = ((a.^3).*((1+7.*(a.^2)).^(-2))).*integral(@(z0)integral_z(z0, a, y, P, A, omega0, to), z0_min, z0_max);

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 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in F_BM_Analytic_Cooling/Integrand (line 56)
    Integrand = ((a.^3).*((1+7.*(a.^2)).^(-2))).*integral(@(z0)integral_z(z0, a, y, P, A, omega0, to), z0_min, z0_max);

Error in F_BM_Analytic_Cooling>@(a)Integrand(a,y,P,A,omega0,to) (line 60)
    IntegrandY = y.^3.*integral(@(a)Integrand(a, y, P, A, omega0, to), 0, 1);

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 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in F_BM_Analytic_Cooling/IntegrandY (line 60)
    IntegrandY = y.^3.*integral(@(a)Integrand(a, y, P, A, omega0, to), 0, 1);

Error in F_BM_Analytic_Cooling>@(y)IntegrandY(y,P,A,omega0,to) (line 27)
    int2 = integral(@(y)IntegrandY(y, P, A, omega0, to),0,1);

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 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in F_BM_Analytic_Cooling (line 27)
    int2 = integral(@(y)IntegrandY(y, P, A, omega0, to),0,1);
    

In order to find out more, I paused the code each time one of the functions is called. That showed me that the first time integral_z is called everything runs smoothly, and y, a and z0 all have the same size (1X150). Then, for a reason I don't understand, integral_z is called a second time with the same y, a only with a different z0, of a different size, (1X90). I don't understand why would integral_z be called again, and I'm really clueless as of what to do.

Any help would be greatly appriciated!

Answer 1

[Walter Roberson]

Each call to integral supplies a vector of values to the function to be called. The vector of values is inherently of inconsistent length. The function must return back a result the same size as the input vector.

You are taking the supplied value, and calling integral() of a different function. But the integral() of the different function supplies its own vector of values of inconsistent length. The sizes of the two vectors disagree and you have problems.

There are three main approaches to take:

1. Code it all using integral3()
2. Code the two outer integrals using the 'arrayvalued', true option
3. Code the two outer integrals to use arrayfun() to apply the functions to one element of the input vector at a time

[Noya Linder]

thank you so much, integral3 does solve the problem! ArrayValued = true took too much time to run, I should've mentioned I tried it previously

Answer 2

[Torsten]

Please supply an executable code in order to test your code and reproduce the error.

Usually, adding 'ArrayValued',true as an option in all your calls to "integral" solves this kind of problem, but it will make your code less efficient.

Consider also a call to "integral3" instead of three calls to "integral".

[Noya Linder]

thank you so much, integral3 does solve the problem! ArrayValued = true took too much time to run, I should've mentioned I tried it previously