Matlab integration error: Integrand output size does not match the input size

I am attempting a 4th order nested integration in Matlab. The function is as follows:

    second_integral = @(r_dash, phi_dash, r, phi) ( (exp(-1i * k * r_dash * cos(phi_dash) * sin(theta)) .* exp(-1i * k * sqrt(r.^2 + r_dash.^2 - 2*( r * r_dash .* cos (phi_dash - phi))))) / (2 .* pi .* sqrt(r.^2 + r_dash.^2 - 2.*( r .* r_dash .* cos (phi_dash - phi))) ) );

With the constants in this particular example being:

    a_1 = 7e-5;
    a_2 = 5e-5;
    k = 1;
    theta = 0; % can also be 0.5, 1, 1.5

The integration is being done as

    integrated_second_integral =integral(@(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi)),0,a_1,'ArrayValued',true);

When the integration command is run I get the following errors:

    Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  2.725290e-19. 
    > In @(r_dash,phi_dash,r,phi)((exp(-1i*k*r_dash*cos(phi_dash)*sin(theta)).*exp(-1i*k*sqrt(r.^2+r_dash.^2-2*(r*r_dash.*cos(phi_dash-phi)))))/(2.*pi.*sqrt(r.^2+r_dash.^2-2.*(r.*r_dash.*cos(phi_dash-phi)))))
      In @(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi)
      In integral3>@(y,z)FUN(x(1)*ones(size(z)),y,z) at 139
      In funfun/private/integral2Calc>integral2t/tensor at 229
      In funfun/private/integral2Calc>integral2t at 56
      In funfun/private/integral2Calc at 10
      In integral3>innerintegral at 138
      In funfun/private/integralCalc>iterateScalarValued at 314
      In funfun/private/integralCalc>vadapt at 133
      In funfun/private/integralCalc at 76
      In integral3 at 122
      In @(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi))
      In funfun/private/integralCalc>iterateArrayValued at 157
      In funfun/private/integralCalc>vadapt at 131
      In funfun/private/integralCalc at 76
      In integral at 89
      In this_matlab_script at 61 
    Error using integral2Calc>integral2t/tensor (line 242)
    Integrand output size does not match the input size.

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

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

    Error in integral3/innerintegral (line 138)
            Q1 = integral2Calc( ...

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

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

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

    Error in integral3 (line 122)
        Q = integralCalc(@innerintegral,xmin,xmax,integralOptions);

    Error in
    @(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi))

    Error in integralCalc/iterateArrayValued (line 157)
                    fxj = FUN(t(1)).*w(1);

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

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

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

    Error in this_matlab_script (line 61)
            integrated_second_integral
            =integral(@(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi)),0,a_1,'ArrayValued',true);

The source of this error is a complete mystery to me. Could someone kindly shed some light on what may be happening here and what I can do to remedy it?

Answer 1

[Mike Hosea]

When you see that error, it means that you used a / instead of ./ . The integrator passed in an array for "vectorized" element-wise evaluation, and instead of just processing it elementwise, your integrand ended trying to solve a linear system with it. This problem is a little bit on the edge of what the routines can solve. I had some problems using the default tolerances, so I had to loosen them a bit and then see how far I could tighten them back up.. After fixing the integrand to be vectorized in all variables, I got this

a_1 = 7e-5;
a_2 = 5e-5;
k = 1;
theta = 0;%can also be 0.5,1,1.5
second_integral=@(r_dash,phi_dash,r,phi)((exp(-1i*k*r_dash.*cos(phi_dash)*sin(theta)).*exp(-1i*k*sqrt(r.^2 + r_dash.^2 - 2*(r.*r_dash.*cos(phi_dash - phi)))))./(2.*pi.*sqrt(r.^2 + r_dash.^2 - 2.*(r.*r_dash.*cos(phi_dash - phi)))));
integrated_second_integral = integralN(second_integral,0,a_1,0,2*pi,0,a_2,0,2*pi,'AbsTol',1e-8,'RelTol',1e-3)

Although theoretically the integral of integral3 approach should work, it's somewhat more efficient to do integral2 of integral2. Setting that up is a little more work because integral2 doesn't have an 'ArrayValued' option, so I used the integralN routine I put on the file exchange which does exactly that.

http://www.mathworks.com/matlabcentral/fileexchange/47919-integraln-m

_

>> integrated_second_integral = integralN(second_integral,0,a_1,0,2*pi,0,a_2,0,2*pi,'AbsTol',1e-8,'RelTol',1e-3)
integrated_second_integral =
      8.479042788173304e-04 - 2.199114856608778e-08i_

[JJ]

Mike, thank you for your answer - I have copied the code and module and it works as expected. I now just need to dig into integralN to understand what it's doing :p Thanks again!